Fix bug in m*L/p calculation introduced by one of the latest commits.
[u/mrichter/AliRoot.git] / JETAN / DEV / AliAnalysisTaskJetCluster.cxx
CommitLineData
d89b8229 1// **************************************
2// Task used for the correction of determiantion of reconstructed jet spectra
3// Compares input (gen) and output (rec) jets
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>
006b2a30 35#include <TF1.h>
d89b8229 36#include <TList.h>
37#include <TLorentzVector.h>
38#include <TClonesArray.h>
39#include "TDatabasePDG.h"
40
41#include "AliAnalysisTaskJetCluster.h"
42#include "AliAnalysisManager.h"
43#include "AliJetFinder.h"
44#include "AliJetHeader.h"
45#include "AliJetReader.h"
46#include "AliESDEvent.h"
47#include "AliAODEvent.h"
48#include "AliAODHandler.h"
49#include "AliAODExtension.h"
50#include "AliAODTrack.h"
51#include "AliAODJet.h"
52#include "AliAODMCParticle.h"
53#include "AliMCEventHandler.h"
54#include "AliMCEvent.h"
55#include "AliStack.h"
56#include "AliGenPythiaEventHeader.h"
57#include "AliJetKineReaderHeader.h"
58#include "AliGenCocktailEventHeader.h"
59#include "AliInputEventHandler.h"
60#include "AliAODJetEventBackground.h"
61
62#include "fastjet/PseudoJet.hh"
63#include "fastjet/ClusterSequenceArea.hh"
64#include "fastjet/AreaDefinition.hh"
65#include "fastjet/JetDefinition.hh"
66// get info on how fastjet was configured
67#include "fastjet/config.h"
68
393a3bd4 69using std::vector;
95808c2b 70
d89b8229 71ClassImp(AliAnalysisTaskJetCluster)
72
73AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
f0659f11 74 //
75 // Destructor
76 //
006b2a30 77
d89b8229 78 delete fRef;
79 delete fRandom;
80
81 if(fTCAJetsOut)fTCAJetsOut->Delete();
82 delete fTCAJetsOut;
006b2a30 83
d89b8229 84 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
85 delete fTCAJetsOutRan;
006b2a30 86
d89b8229 87 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
88 delete fTCARandomConesOut;
006b2a30 89
d89b8229 90 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
91 delete fTCARandomConesOutRan;
006b2a30 92
d89b8229 93 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
94 delete fAODJetBackgroundOut;
95}
96
97AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster():
98 AliAnalysisTaskSE(),
99 fAOD(0x0),
100 fAODExtension(0x0),
101 fRef(new TRefArray),
102 fUseAODTrackInput(kFALSE),
103 fUseAODMCInput(kFALSE),
104 fUseBackgroundCalc(kFALSE),
105 fEventSelection(kFALSE),
106 fFilterMask(0),
d7396b04 107 fFilterMaskBestPt(0),
d89b8229 108 fFilterType(0),
109 fJetTypes(kJet),
110 fTrackTypeRec(kTrackUndef),
111 fTrackTypeGen(kTrackUndef),
112 fNSkipLeadingRan(0),
113 fNSkipLeadingCone(0),
114 fNRandomCones(0),
115 fAvgTrials(1),
116 fExternalWeight(1),
117 fTrackEtaWindow(0.9),
118 fRecEtaWindow(0.5),
119 fTrackPtCut(0.),
120 fJetOutputMinPt(0.150),
121 fMaxTrackPtInJet(100.),
122 fJetTriggerPtCut(0),
123 fVtxZCut(8),
124 fVtxR2Cut(1),
125 fCentCutUp(0),
126 fCentCutLo(0),
127 fNonStdBranch(""),
128 fBackgroundBranch(""),
129 fNonStdFile(""),
006b2a30 130 fMomResH1(0x0),
131 fMomResH2(0x0),
132 fMomResH3(0x0),
133 fMomResH1Fit(0x0),
134 fMomResH2Fit(0x0),
135 fMomResH3Fit(0x0),
136 fhEffH1(0x0),
137 fhEffH2(0x0),
138 fhEffH3(0x0),
fa7c34ba 139 fUseTrPtResolutionSmearing(kFALSE),
006b2a30 140 fUseDiceEfficiency(kFALSE),
4bd66a92 141 fDiceEfficiencyMinPt(0.),
fa7c34ba 142 fUseTrPtResolutionFromOADB(kFALSE),
143 fUseTrEfficiencyFromOADB(kFALSE),
144 fPathTrPtResolution(""),
145 fPathTrEfficiency(""),
146 fChangeEfficiencyFraction(0.),
d89b8229 147 fRparam(1.0),
148 fAlgorithm(fastjet::kt_algorithm),
149 fStrategy(fastjet::Best),
150 fRecombScheme(fastjet::BIpt_scheme),
151 fAreaType(fastjet::active_area),
152 fGhostArea(0.01),
153 fActiveAreaRepeats(1),
154 fGhostEtamax(1.5),
155 fTCAJetsOut(0x0),
156 fTCAJetsOutRan(0x0),
157 fTCARandomConesOut(0x0),
158 fTCARandomConesOutRan(0x0),
159 fAODJetBackgroundOut(0x0),
160 fRandom(0),
161 fh1Xsec(0x0),
162 fh1Trials(0x0),
163 fh1PtHard(0x0),
164 fh1PtHardNoW(0x0),
165 fh1PtHardTrials(0x0),
166 fh1NJetsRec(0x0),
167 fh1NConstRec(0x0),
168 fh1NConstLeadingRec(0x0),
169 fh1PtJetsRecIn(0x0),
170 fh1PtJetsLeadingRecIn(0x0),
171 fh1PtJetConstRec(0x0),
172 fh1PtJetConstLeadingRec(0x0),
173 fh1PtTracksRecIn(0x0),
174 fh1PtTracksLeadingRecIn(0x0),
175 fh1NJetsRecRan(0x0),
176 fh1NConstRecRan(0x0),
177 fh1PtJetsLeadingRecInRan(0x0),
178 fh1NConstLeadingRecRan(0x0),
179 fh1PtJetsRecInRan(0x0),
180 fh1PtTracksGenIn(0x0),
181 fh1Nch(0x0),
182 fh1CentralityPhySel(0x0),
183 fh1Centrality(0x0),
184 fh1CentralitySelect(0x0),
185 fh1ZPhySel(0x0),
186 fh1Z(0x0),
187 fh1ZSelect(0x0),
188 fh2NRecJetsPt(0x0),
189 fh2NRecTracksPt(0x0),
190 fh2NConstPt(0x0),
191 fh2NConstLeadingPt(0x0),
192 fh2JetPhiEta(0x0),
193 fh2LeadingJetPhiEta(0x0),
194 fh2JetEtaPt(0x0),
195 fh2LeadingJetEtaPt(0x0),
196 fh2TrackEtaPt(0x0),
197 fh2LeadingTrackEtaPt(0x0),
198 fh2JetsLeadingPhiEta(0x0),
199 fh2JetsLeadingPhiPt(0x0),
200 fh2TracksLeadingPhiEta(0x0),
201 fh2TracksLeadingPhiPt(0x0),
202 fh2TracksLeadingJetPhiPt(0x0),
203 fh2JetsLeadingPhiPtW(0x0),
204 fh2TracksLeadingPhiPtW(0x0),
205 fh2TracksLeadingJetPhiPtW(0x0),
206 fh2NRecJetsPtRan(0x0),
207 fh2NConstPtRan(0x0),
208 fh2NConstLeadingPtRan(0x0),
209 fh2PtNch(0x0),
210 fh2PtNchRan(0x0),
211 fh2PtNchN(0x0),
212 fh2PtNchNRan(0x0),
213 fh2TracksLeadingJetPhiPtRan(0x0),
214 fh2TracksLeadingJetPhiPtWRan(0x0),
c0fdfe53 215 fh3CentvsRhoLeadingTrackCorr(0x0),
216 fh3CentvsSigmaLeadingTrackCorr(0x0),
217 fh3MultvsRhoLeadingTrackCorr(0x0),
218 fh3MultvsSigmaLeadingTrackCorr(0x0),
006b2a30 219 fh2PtGenPtSmeared(0x0),
220 fp1Efficiency(0x0),
221 fp1PtResolution(0x0),
d89b8229 222 fHistList(0x0)
223{
f0659f11 224 //
225 // Constructor
226 //
227
d89b8229 228 for(int i = 0;i<3;i++){
229 fh1BiARandomCones[i] = 0;
230 fh1BiARandomConesRan[i] = 0;
231 }
232 for(int i = 0;i<kMaxCent;i++){
233 fh2JetsLeadingPhiPtC[i] = 0;
234 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
235 fh2TracksLeadingJetPhiPtC[i] = 0;
236 fh2TracksLeadingJetPhiPtWC[i] = 0;
237 }
238}
239
240AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
241 AliAnalysisTaskSE(name),
242 fAOD(0x0),
243 fAODExtension(0x0),
244 fRef(new TRefArray),
245 fUseAODTrackInput(kFALSE),
246 fUseAODMCInput(kFALSE),
247 fUseBackgroundCalc(kFALSE),
d7396b04 248 fEventSelection(kFALSE), fFilterMask(0),
249 fFilterMaskBestPt(0),
d89b8229 250 fFilterType(0),
251 fJetTypes(kJet),
252 fTrackTypeRec(kTrackUndef),
253 fTrackTypeGen(kTrackUndef),
254 fNSkipLeadingRan(0),
255 fNSkipLeadingCone(0),
256 fNRandomCones(0),
257 fAvgTrials(1),
258 fExternalWeight(1),
259 fTrackEtaWindow(0.9),
260 fRecEtaWindow(0.5),
261 fTrackPtCut(0.),
262 fJetOutputMinPt(0.150),
263 fMaxTrackPtInJet(100.),
264 fJetTriggerPtCut(0),
265 fVtxZCut(8),
266 fVtxR2Cut(1),
267 fCentCutUp(0),
268 fCentCutLo(0),
269 fNonStdBranch(""),
270 fBackgroundBranch(""),
271 fNonStdFile(""),
006b2a30 272 fMomResH1(0x0),
273 fMomResH2(0x0),
274 fMomResH3(0x0),
275 fMomResH1Fit(0x0),
276 fMomResH2Fit(0x0),
277 fMomResH3Fit(0x0),
278 fhEffH1(0x0),
279 fhEffH2(0x0),
280 fhEffH3(0x0),
fa7c34ba 281 fUseTrPtResolutionSmearing(kFALSE),
006b2a30 282 fUseDiceEfficiency(kFALSE),
4bd66a92 283 fDiceEfficiencyMinPt(0.),
fa7c34ba 284 fUseTrPtResolutionFromOADB(kFALSE),
285 fUseTrEfficiencyFromOADB(kFALSE),
286 fPathTrPtResolution(""),
287 fPathTrEfficiency(""),
288 fChangeEfficiencyFraction(0.),
d89b8229 289 fRparam(1.0),
290 fAlgorithm(fastjet::kt_algorithm),
291 fStrategy(fastjet::Best),
292 fRecombScheme(fastjet::BIpt_scheme),
293 fAreaType(fastjet::active_area),
294 fGhostArea(0.01),
295 fActiveAreaRepeats(1),
296 fGhostEtamax(1.5),
297 fTCAJetsOut(0x0),
298 fTCAJetsOutRan(0x0),
299 fTCARandomConesOut(0x0),
300 fTCARandomConesOutRan(0x0),
301 fAODJetBackgroundOut(0x0),
302 fRandom(0),
303 fh1Xsec(0x0),
304 fh1Trials(0x0),
305 fh1PtHard(0x0),
306 fh1PtHardNoW(0x0),
307 fh1PtHardTrials(0x0),
308 fh1NJetsRec(0x0),
309 fh1NConstRec(0x0),
310 fh1NConstLeadingRec(0x0),
311 fh1PtJetsRecIn(0x0),
312 fh1PtJetsLeadingRecIn(0x0),
313 fh1PtJetConstRec(0x0),
314 fh1PtJetConstLeadingRec(0x0),
315 fh1PtTracksRecIn(0x0),
316 fh1PtTracksLeadingRecIn(0x0),
317 fh1NJetsRecRan(0x0),
318 fh1NConstRecRan(0x0),
319 fh1PtJetsLeadingRecInRan(0x0),
320 fh1NConstLeadingRecRan(0x0),
321 fh1PtJetsRecInRan(0x0),
322 fh1PtTracksGenIn(0x0),
323 fh1Nch(0x0),
324 fh1CentralityPhySel(0x0),
325 fh1Centrality(0x0),
326 fh1CentralitySelect(0x0),
327 fh1ZPhySel(0x0),
328 fh1Z(0x0),
329 fh1ZSelect(0x0),
330 fh2NRecJetsPt(0x0),
331 fh2NRecTracksPt(0x0),
332 fh2NConstPt(0x0),
333 fh2NConstLeadingPt(0x0),
334 fh2JetPhiEta(0x0),
335 fh2LeadingJetPhiEta(0x0),
336 fh2JetEtaPt(0x0),
337 fh2LeadingJetEtaPt(0x0),
338 fh2TrackEtaPt(0x0),
339 fh2LeadingTrackEtaPt(0x0),
340 fh2JetsLeadingPhiEta(0x0),
341 fh2JetsLeadingPhiPt(0x0),
342 fh2TracksLeadingPhiEta(0x0),
343 fh2TracksLeadingPhiPt(0x0),
344 fh2TracksLeadingJetPhiPt(0x0),
345 fh2JetsLeadingPhiPtW(0x0),
346 fh2TracksLeadingPhiPtW(0x0),
347 fh2TracksLeadingJetPhiPtW(0x0),
348 fh2NRecJetsPtRan(0x0),
349 fh2NConstPtRan(0x0),
350 fh2NConstLeadingPtRan(0x0),
351 fh2PtNch(0x0),
352 fh2PtNchRan(0x0),
353 fh2PtNchN(0x0),
354 fh2PtNchNRan(0x0),
355 fh2TracksLeadingJetPhiPtRan(0x0),
356 fh2TracksLeadingJetPhiPtWRan(0x0),
c0fdfe53 357 fh3CentvsRhoLeadingTrackCorr(0x0),
358 fh3CentvsSigmaLeadingTrackCorr(0x0),
359 fh3MultvsRhoLeadingTrackCorr(0x0),
360 fh3MultvsSigmaLeadingTrackCorr(0x0),
006b2a30 361 fh2PtGenPtSmeared(0x0),
362 fp1Efficiency(0x0),
363 fp1PtResolution(0x0),
d89b8229 364 fHistList(0x0)
365{
f0659f11 366 //
367 // named ctor
368 //
006b2a30 369
d89b8229 370 for(int i = 0;i<3;i++){
371 fh1BiARandomCones[i] = 0;
372 fh1BiARandomConesRan[i] = 0;
373 }
374 for(int i = 0;i<kMaxCent;i++){
375 fh2JetsLeadingPhiPtC[i] = 0;
376 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
377 fh2TracksLeadingJetPhiPtC[i] = 0;
378 fh2TracksLeadingJetPhiPtWC[i] = 0;
379 }
380 DefineOutput(1, TList::Class());
381}
382
383
384
385Bool_t AliAnalysisTaskJetCluster::Notify()
386{
387 //
388 // Implemented Notify() to read the cross sections
389 // and number of trials from pyxsec.root
390 //
391 return kTRUE;
392}
393
394void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
395{
396
397 //
398 // Create the output container
399 //
400
401 fRandom = new TRandom3(0);
402
403
404 // Connect the AOD
405
406
407 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
408
409
410
411 if(fNonStdBranch.Length()!=0)
412 {
413 // only create the output branch if we have a name
414 // Create a new branch for jets...
415 // -> cleared in the UserExec....
416 // here we can also have the case that the brnaches are written to a separate file
417
418 if(fJetTypes&kJet){
419 fTCAJetsOut = new TClonesArray("AliAODJet", 0);
420 fTCAJetsOut->SetName(fNonStdBranch.Data());
421 AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
422 }
006b2a30 423
d89b8229 424 if(fJetTypes&kJetRan){
425 fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
426 fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
fa7c34ba 427 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing)
428 fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dFr%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
d89b8229 429 AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
430 }
431
432 if(fUseBackgroundCalc){
433 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
434 fAODJetBackgroundOut = new AliAODJetEventBackground();
435 fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
fa7c34ba 436 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing)
437 fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dFr%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
006b2a30 438
d89b8229 439 AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
440 }
441 }
442 // create the branch for the random cones with the same R
443 TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
fa7c34ba 444 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing)
445 cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone);
006b2a30 446
d89b8229 447 if(fNRandomCones>0){
448 if(fJetTypes&kRC){
449 if(!AODEvent()->FindListObject(cName.Data())){
450 fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
451 fTCARandomConesOut->SetName(cName.Data());
452 AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
453 }
454 }
455 // create the branch with the random for the random cones on the random event
456 if(fJetTypes&kRCRan){
457 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
458 if(!AODEvent()->FindListObject(cName.Data())){
459 fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
460 fTCARandomConesOutRan->SetName(cName.Data());
461 AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
462 }
463 }
464 }
465
466 if(fNonStdFile.Length()!=0){
467 //
468 // case that we have an AOD extension we need to fetch the jets from the extended output
469 // we identify the extension aod event by looking for the branchname
470 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
471 // case that we have an AOD extension we need can fetch the background maybe from the extended output
472 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
473 }
474 }
475
476
477 if(!fHistList)fHistList = new TList();
478 fHistList->SetOwner();
479 PostData(1, fHistList); // post data in any case once
480
481 Bool_t oldStatus = TH1::AddDirectoryStatus();
482 TH1::AddDirectory(kFALSE);
483
484 //
485 // Histogram
486
487 const Int_t nBinPt = 100;
488 Double_t binLimitsPt[nBinPt+1];
489 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
490 if(iPt == 0){
491 binLimitsPt[iPt] = 0.0;
492 }
493 else {// 1.0
494 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
495 }
496 }
497
498 const Int_t nBinPhi = 90;
499 Double_t binLimitsPhi[nBinPhi+1];
500 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
501 if(iPhi==0){
502 binLimitsPhi[iPhi] = -1.*TMath::Pi();
503 }
504 else{
505 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
506 }
507 }
508
509
510
511 const Int_t nBinEta = 40;
512 Double_t binLimitsEta[nBinEta+1];
513 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
514 if(iEta==0){
515 binLimitsEta[iEta] = -2.0;
516 }
517 else{
518 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
519 }
520 }
521
522 const int nChMax = 5000;
523
524 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
525 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
526
527 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
528 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
529
530
531 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
532 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
533
534 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
535 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
536 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
537 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
538
539
540 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
541 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
542 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
543
544 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
545 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
546 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
547 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
548 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
549 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
550 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
551 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
552 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
553 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
554
555 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
556 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
557 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
558
559 fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
560 fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
561 fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
562
563 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
564 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
565 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
566 //
567
568
569 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
570 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
571 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
572 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
573
574 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
575 fh2PtNchRan = new TH2F("fh2PtNchRan","p_T of cluster vs. multiplicity ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
576 fh2PtNchN = new TH2F("fh2PtNchN","p_T of cluster vs. multiplicity N weighted; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
577 fh2PtNchNRan = new TH2F("fh2PtNchNRan","p_T of cluster vs. multiplicity N weighted ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
578
579
580
581 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
582 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
583 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
584 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
585
586 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
587 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
588 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
589 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
590
591 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
592 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
593 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
594 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
595
596
597
598 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
599 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
600 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
601 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
602 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
603 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
604 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
605 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
606 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
607 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
608 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
609 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
610
611 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
612 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
613 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
614 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
615
616 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
617 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
618 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
619 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
620
c0fdfe53 621 fh3CentvsRhoLeadingTrackCorr = new TH3F("fh3CentvsRhoLeadingTrackCorr","centrality vs background density; centrality; #rho; Quadrant", 100,0.,100.,600,0.,300.,4,0.,4.);
622 fh3CentvsSigmaLeadingTrackCorr = new TH3F("fh3CentvsSigmaLeadingTrackCorr","centrality vs backgroun sigma; centrality; #sigma(#rho); Quadrant",100,0.,100.,500,0.,50.,4,0.,4.);
623 fHistList->Add(fh3CentvsRhoLeadingTrackCorr);
624 fHistList->Add(fh3CentvsSigmaLeadingTrackCorr);
625
626 fh3MultvsRhoLeadingTrackCorr = new TH3F("fh3MultvsRhoLeadingTrackCorr","mult vs background density; N_{tracks}; #rho; Quadrant", 100,0.,5000.,600,0.,300.,4,0.,4.);
627 fh3MultvsSigmaLeadingTrackCorr = new TH3F("fh3MultvsSigmaLeadingTrackCorr","mult vs background sigma; N_{tracks}; #sigma(#rho); Quadrant",100,0.,5000.,500,0.,50.,4,0.,4.);
628 fHistList->Add(fh3MultvsRhoLeadingTrackCorr);
629 fHistList->Add(fh3MultvsSigmaLeadingTrackCorr);
630
006b2a30 631 //Detector level effects histos
632 fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
633
634 fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
635 fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
636
637 fHistList->Add(fh2PtGenPtSmeared);
638 fHistList->Add(fp1Efficiency);
639 fHistList->Add(fp1PtResolution);
d89b8229 640
641 if(fNRandomCones>0&&fUseBackgroundCalc){
642 for(int i = 0;i<3;i++){
643 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
644 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
645 }
646 }
647
648 for(int i = 0;i < kMaxCent;i++){
649 fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
650 fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
651 fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
652 fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
653 }
654
655 const Int_t saveLevel = 3; // large save level more histos
656 if(saveLevel>0){
657 fHistList->Add(fh1Xsec);
658 fHistList->Add(fh1Trials);
659
660 fHistList->Add(fh1NJetsRec);
661 fHistList->Add(fh1NConstRec);
662 fHistList->Add(fh1NConstLeadingRec);
663 fHistList->Add(fh1PtJetsRecIn);
664 fHistList->Add(fh1PtJetsLeadingRecIn);
665 fHistList->Add(fh1PtTracksRecIn);
666 fHistList->Add(fh1PtTracksLeadingRecIn);
667 fHistList->Add(fh1PtJetConstRec);
668 fHistList->Add(fh1PtJetConstLeadingRec);
669 fHistList->Add(fh1NJetsRecRan);
670 fHistList->Add(fh1NConstRecRan);
671 fHistList->Add(fh1PtJetsLeadingRecInRan);
672 fHistList->Add(fh1NConstLeadingRecRan);
673 fHistList->Add(fh1PtJetsRecInRan);
674 fHistList->Add(fh1Nch);
675 fHistList->Add(fh1Centrality);
676 fHistList->Add(fh1CentralitySelect);
677 fHistList->Add(fh1CentralityPhySel);
678 fHistList->Add(fh1Z);
679 fHistList->Add(fh1ZSelect);
680 fHistList->Add(fh1ZPhySel);
681 if(fNRandomCones>0&&fUseBackgroundCalc){
682 for(int i = 0;i<3;i++){
683 fHistList->Add(fh1BiARandomCones[i]);
684 fHistList->Add(fh1BiARandomConesRan[i]);
685 }
686 }
687 for(int i = 0;i < kMaxCent;i++){
688 fHistList->Add(fh2JetsLeadingPhiPtC[i]);
689 fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
690 fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
691 fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
692 }
693
694 fHistList->Add(fh2NRecJetsPt);
695 fHistList->Add(fh2NRecTracksPt);
696 fHistList->Add(fh2NConstPt);
697 fHistList->Add(fh2NConstLeadingPt);
698 fHistList->Add(fh2PtNch);
699 fHistList->Add(fh2PtNchRan);
700 fHistList->Add(fh2PtNchN);
701 fHistList->Add(fh2PtNchNRan);
702 fHistList->Add(fh2JetPhiEta);
703 fHistList->Add(fh2LeadingJetPhiEta);
704 fHistList->Add(fh2JetEtaPt);
705 fHistList->Add(fh2LeadingJetEtaPt);
706 fHistList->Add(fh2TrackEtaPt);
707 fHistList->Add(fh2LeadingTrackEtaPt);
708 fHistList->Add(fh2JetsLeadingPhiEta );
709 fHistList->Add(fh2JetsLeadingPhiPt);
710 fHistList->Add(fh2TracksLeadingPhiEta);
711 fHistList->Add(fh2TracksLeadingPhiPt);
712 fHistList->Add(fh2TracksLeadingJetPhiPt);
713 fHistList->Add(fh2JetsLeadingPhiPtW);
714 fHistList->Add(fh2TracksLeadingPhiPtW);
715 fHistList->Add(fh2TracksLeadingJetPhiPtW);
716 fHistList->Add(fh2NRecJetsPtRan);
717 fHistList->Add(fh2NConstPtRan);
718 fHistList->Add(fh2NConstLeadingPtRan);
719 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
720 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
721 }
722
723 // =========== Switch on Sumw2 for all histos ===========
724 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
725 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
726 if (h1){
727 h1->Sumw2();
728 continue;
729 }
730 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
731 if(hn)hn->Sumw2();
732 }
733 TH1::AddDirectory(oldStatus);
734}
735
fa7c34ba 736void AliAnalysisTaskJetCluster::LocalInit()
d89b8229 737{
738 //
739 // Initialization
740 //
741
742 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
743
fa7c34ba 744 if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB();
745 if(fUseTrEfficiencyFromOADB) LoadTrEfficiencyRootFileFromOADB();
746
747
006b2a30 748 FitMomentumResolution();
749
d89b8229 750}
751
752void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
753{
754
755 // handle and reset the output jet branch
756
757 if(fTCAJetsOut)fTCAJetsOut->Delete();
758 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
759 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
760 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
761 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
762
763 AliAODJetEventBackground* externalBackground = 0;
764 if(!externalBackground&&fBackgroundBranch.Length()){
765 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
766 if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
767 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
768 }
769 //
770 // Execute analysis for current event
771 //
772 AliESDEvent *fESD = 0;
773 if(fUseAODTrackInput){
774 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
775 if(!fAOD){
776 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
777 return;
778 }
006b2a30 779 // fetch the header
d89b8229 780 }
781 else{
782 // assume that the AOD is in the general output...
783 fAOD = AODEvent();
784 if(!fAOD){
785 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
786 return;
787 }
788 if(fDebug>0){
789 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
790 }
791 }
006b2a30 792
793 //Check if information is provided detector level effects
fa7c34ba 794 if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE;
006b2a30 795 if(!fhEffH1 || !fhEffH2 || !fhEffH3) fUseDiceEfficiency = kFALSE;
d89b8229 796
797 Bool_t selectEvent = false;
798 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
799
800 Float_t cent = 0;
801 Float_t zVtx = 0;
802 Int_t cenClass = -1;
803 if(fAOD){
804 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
805 TString vtxTitle(vtxAOD->GetTitle());
806 zVtx = vtxAOD->GetZ();
807
808 cent = fAOD->GetHeader()->GetCentrality();
809 if(cent<10)cenClass = 0;
810 else if(cent<30)cenClass = 1;
811 else if(cent<50)cenClass = 2;
812 else if(cent<80)cenClass = 3;
813 if(physicsSelection){
814 fh1CentralityPhySel->Fill(cent);
815 fh1ZPhySel->Fill(zVtx);
816 }
817
818 if(fEventSelection){
819 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
820 Float_t yvtx = vtxAOD->GetY();
821 Float_t xvtx = vtxAOD->GetX();
822 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
823 if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
824 if(physicsSelection){
825 selectEvent = true;
826 }
827 }
828 }
829 if(fCentCutUp>0){
830 if(cent<fCentCutLo||cent>fCentCutUp){
831 selectEvent = false;
832 }
833 }
834 }else{
835 selectEvent = true;
836 }
837 }
838
839
840 if(!selectEvent){
841 PostData(1, fHistList);
842 return;
843 }
844 fh1Centrality->Fill(cent);
845 fh1Z->Fill(zVtx);
846 fh1Trials->Fill("#sum{ntrials}",1);
847
848
849 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
850
851 // ==== General variables needed
852
853
854
855 // we simply fetch the tracks/mc particles as a list of AliVParticles
856
857 TList recParticles;
858 TList genParticles;
859
860 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
861 Float_t nCh = recParticles.GetEntries();
862 fh1Nch->Fill(nCh);
863 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
864 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
865 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
866
867 // find the jets....
868
869 vector<fastjet::PseudoJet> inputParticlesRec;
870 vector<fastjet::PseudoJet> inputParticlesRecRan;
871
872 // Generate the random cones
873
874 AliAODJet vTmpRan(1,0,0,1);
875 for(int i = 0; i < recParticles.GetEntries(); i++){
876 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
006b2a30 877
d89b8229 878 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
879 // we take total momentum here
006b2a30 880
fa7c34ba 881 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
006b2a30 882 //Add particles to fastjet in case we are not running toy model
883 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
884 jInp.set_user_index(i);
885 inputParticlesRec.push_back(jInp);
886 }
887 else if(fUseDiceEfficiency) {
888
889 // Dice to decide if particle is kept or not - toy model for efficiency
890 //
891 Double_t rnd = fRandom->Uniform(1.);
892 Double_t pT = vp->Pt();
893 Double_t eff[3] = {0.};
894 Double_t pTtmp = pT;
895 if(pT>10.) pTtmp = 10.;
896 Double_t eff1 = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
897 Double_t eff2 = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
898 Double_t eff3 = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));
899 Int_t cat[3] = {0};
900 //Sort efficiencies from large to small
901 if(eff1>eff2 && eff1>eff3) {
902 eff[0] = eff1;
903 cat[0] = 1;
904 if(eff2>eff3) {
905 eff[1] = eff2;
906 eff[2] = eff3;
907 cat[1] = 2;
908 cat[2] = 3;
909 } else {
910 eff[1] = eff3;
911 eff[2] = eff2;
912 cat[1] = 3;
913 cat[2] = 2;
914 }
915 }
916 else if(eff2>eff1 && eff2>eff3) {
917 eff[0] = eff2;
918 cat[0] = 2;
919 if(eff1>eff3) {
920 eff[1] = eff1;
921 eff[2] = eff3;
922 cat[1] = 1;
923 cat[2] = 3;
924 } else {
925 eff[1] = eff3;
926 eff[2] = eff1;
927 cat[1] = 3;
928 cat[2] = 1;
929 }
930 }
931 else if(eff3>eff1 && eff3>eff2) {
932 eff[0] = eff3;
933 cat[0] = 3;
934 if(eff1>eff2) {
935 eff[1] = eff1;
936 eff[2] = eff2;
937 cat[1] = 1;
938 cat[2] = 2;
939 } else {
940 eff[1] = eff2;
941 eff[2] = eff1;
942 cat[1] = 2;
943 cat[2] = 1;
944 }
945 }
946
947 Double_t sumEff = eff[0]+eff[1]+eff[2];
948 fp1Efficiency->Fill(vp->Pt(),sumEff);
4bd66a92 949 if(rnd>sumEff && pT > fDiceEfficiencyMinPt) continue;
006b2a30 950
fa7c34ba 951 if(fUseTrPtResolutionSmearing) {
006b2a30 952 //Smear momentum of generated particle
953 Double_t smear = 1.;
954 //Select hybrid track category
955 if(rnd<=eff[2])
956 smear = GetMomentumSmearing(cat[2],pT);
957 else if(rnd<=(eff[2]+eff[1]))
958 smear = GetMomentumSmearing(cat[1],pT);
959 else
960 smear = GetMomentumSmearing(cat[0],pT);
961
962 fp1PtResolution->Fill(vp->Pt(),smear);
963
964 Double_t sigma = vp->Pt()*smear;
965 Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
966
967 Double_t phi = vp->Phi();
968 Double_t theta = 2.*TMath::ATan(TMath::Exp(-1.*(vp->Eta())));
969 Double_t pX = pTrec * TMath::Cos(phi);
970 Double_t pY = pTrec * TMath::Sin(phi);
971 Double_t pZ = pTrec/TMath::Tan(theta);
972 Double_t p=TMath::Sqrt(pTrec*pTrec+pZ*pZ);
973
974 fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);
975
976 fastjet::PseudoJet jInp(pX,pY,pZ,p);
977 jInp.set_user_index(i);
978 inputParticlesRec.push_back(jInp);
979
980 }
981 else {
982 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
983 jInp.set_user_index(i);
984 inputParticlesRec.push_back(jInp);
985
986 }
987
988 }
d89b8229 989
990 // the randomized input changes eta and phi, but keeps the p_T
991 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
992 Double_t pT = vp->Pt();
993 Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
994 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
995
996 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
997 Double_t pZ = pT/TMath::Tan(theta);
998
999 Double_t pX = pT * TMath::Cos(phi);
1000 Double_t pY = pT * TMath::Sin(phi);
1001 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
1002 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
1003
1004 jInpRan.set_user_index(i);
1005 inputParticlesRecRan.push_back(jInpRan);
1006 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
1007 }
1008
1009 // fill the tref array, only needed when we write out jets
1010 if(fTCAJetsOut){
1011 if(i == 0){
1012 fRef->Delete(); // make sure to delete before placement new...
fa7c34ba 1013 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
006b2a30 1014 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
1015 }
d89b8229 1016 }
fa7c34ba 1017 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ...
d89b8229 1018 }
1019 }// recparticles
1020
1021 if(inputParticlesRec.size()==0){
1022 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
1023 PostData(1, fHistList);
1024 return;
1025 }
1026
1027 // run fast jet
1028 // employ setters for these...
1029
1030
1031 // now create the object that holds info about ghosts
1032 /*
1033 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
1034 // reduce CPU time...
1035 fGhostArea = 0.5;
1036 fActiveAreaRepeats = 0;
1037 }
1038 */
1039
1040 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
1041 fastjet::AreaType areaType = fastjet::active_area;
1042 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
1043 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
1044 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
1045
1046 //range where to compute background
1047 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
1048 phiMin = 0;
1049 phiMax = 2*TMath::Pi();
1050 rapMax = fGhostEtamax - fRparam;
1051 rapMin = - fGhostEtamax + fRparam;
1052 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
1053
1054
1055 const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
1056 const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
1057
1058
1059 fh1NJetsRec->Fill(sortedJets.size());
1060
1061 // loop over all jets an fill information, first on is the leading jet
1062
1063 Int_t nRecOver = inclusiveJets.size();
1064 Int_t nRec = inclusiveJets.size();
1065 if(inclusiveJets.size()>0){
1066 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
1067 Double_t area = clustSeq.area(sortedJets[0]);
1068 leadingJet.SetEffArea(area,0);
1069 Float_t pt = leadingJet.Pt();
1070 Int_t nAodOutJets = 0;
1071 Int_t nAodOutTracks = 0;
1072 AliAODJet *aodOutJet = 0;
1073
1074 Int_t iCount = 0;
1075 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
1076 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
1077 while(pt<ptCut&&iCount<nRec){
1078 nRecOver--;
1079 iCount++;
1080 if(iCount<nRec){
1081 pt = sortedJets[iCount].perp();
1082 }
1083 }
1084 if(nRecOver<=0)break;
1085 fh2NRecJetsPt->Fill(ptCut,nRecOver);
1086 }
1087 Float_t phi = leadingJet.Phi();
1088 if(phi<0)phi+=TMath::Pi()*2.;
1089 Float_t eta = leadingJet.Eta();
1090 Float_t pTback = 0;
1091 if(externalBackground){
1092 // carefull has to be filled in a task before
006b2a30 1093 // todo, ReArrange to the botom
83decb5f 1094 pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
d89b8229 1095 }
1096 pt = leadingJet.Pt() - pTback;
1097 // correlation of leading jet with tracks
1098 TIterator *recIter = recParticles.MakeIterator();
1099 recIter->Reset();
1100 AliVParticle *tmpRecTrack = 0;
1101 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
1102 Float_t tmpPt = tmpRecTrack->Pt();
1103 // correlation
1104 Float_t tmpPhi = tmpRecTrack->Phi();
1105 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1106 Float_t dPhi = phi - tmpPhi;
1107 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1108 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1109 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
1110 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
1111 if(cenClass>=0){
1112 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
1113 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1114 }
1115
1116 }
1117
1118
1119 TLorentzVector vecareab;
1120 for(int j = 0; j < nRec;j++){
1121 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
1122 aodOutJet = 0;
1123 nAodOutTracks = 0;
1124 Float_t tmpPt = tmpRec.Pt();
1125
1126 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
1127 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
1128 aodOutJet->GetRefTracks()->Clear();
1129 Double_t area1 = clustSeq.area(sortedJets[j]);
1130 aodOutJet->SetEffArea(area1,0);
1131 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
1132 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
1133 aodOutJet->SetVectorAreaCharged(&vecareab);
1134 }
1135
1136
1137 Float_t tmpPtBack = 0;
1138 if(externalBackground){
1139 // carefull has to be filled in a task before
1140 // todo, ReArrange to the botom
1141 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
1142 }
1143 tmpPt = tmpPt - tmpPtBack;
1144 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
1145
1146 fh1PtJetsRecIn->Fill(tmpPt);
1147 // Fill Spectra with constituentsemacs
1148 const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
1149
1150 fh1NConstRec->Fill(constituents.size());
1151 fh2PtNch->Fill(nCh,tmpPt);
1152 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
1153 fh2NConstPt->Fill(tmpPt,constituents.size());
1154 // loop over constiutents and fill spectrum
1155
d7396b04 1156 AliVParticle *partLead = 0;
1157 Float_t ptLead = -1;
1158
d89b8229 1159 for(unsigned int ic = 0; ic < constituents.size();ic++){
1160 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
006b2a30 1161 if(!part) continue;
d89b8229 1162 fh1PtJetConstRec->Fill(part->Pt());
1163 if(aodOutJet){
fa7c34ba 1164 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
d7396b04 1165 if(part->Pt()>fMaxTrackPtInJet){
1166 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1167 }
1168 }
1169 if(part->Pt()>ptLead){
95808c2b 1170 ptLead = part->Pt();
d7396b04 1171 partLead = part;
d89b8229 1172 }
1173 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1174 }
d7396b04 1175
1176 AliAODTrack *aodT = 0;
47373258 1177 if(partLead){
1178 if(aodOutJet){
1179 //set pT of leading constituent of jet
1180 aodOutJet->SetPtLeading(partLead->Pt());
1181 aodT = dynamic_cast<AliAODTrack*>(partLead);
1182 if(aodT){
1183 if(aodT->TestFilterBit(fFilterMaskBestPt)){
1184 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
1185 }
d7396b04 1186 }
1187 }
1188 }
47373258 1189
d89b8229 1190 // correlation
1191 Float_t tmpPhi = tmpRec.Phi();
1192 Float_t tmpEta = tmpRec.Eta();
1193 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1194 if(j==0){
1195 fh1PtJetsLeadingRecIn->Fill(tmpPt);
1196 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1197 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1198 fh1NConstLeadingRec->Fill(constituents.size());
1199 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1200 continue;
1201 }
1202 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1203 fh2JetEtaPt->Fill(tmpEta,tmpPt);
1204 Float_t dPhi = phi - tmpPhi;
1205 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1206 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1207 Float_t dEta = eta - tmpRec.Eta();
1208 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1209 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1210 if(cenClass>=0){
1211 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1212 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1213 }
1214 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1215 }// loop over reconstructed jets
1216 delete recIter;
1217
1218
1219
1220 // Add the random cones...
1221 if(fNRandomCones>0&&fTCARandomConesOut){
1222 // create a random jet within the acceptance
1223 Double_t etaMax = fTrackEtaWindow - fRparam;
1224 Int_t nCone = 0;
1225 Int_t nConeRan = 0;
1226 Double_t pTC = 1; // small number
1227 for(int ir = 0;ir < fNRandomCones;ir++){
1228 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1229 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1230 // massless jet
1231 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1232 Double_t pZC = pTC/TMath::Tan(thetaC);
1233 Double_t pXC = pTC * TMath::Cos(phiC);
1234 Double_t pYC = pTC * TMath::Sin(phiC);
1235 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1236 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
1237 bool skip = false;
1238 for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1239 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1240 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1241 skip = true;
1242 break;
1243 }
1244 }
1245 // test for overlap with previous cones to avoid double counting
1246 for(int iic = 0;iic<ir;iic++){
1247 AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1248 if(iicone){
1249 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1250 skip = true;
1251 break;
1252 }
1253 }
1254 }
1255 if(skip)continue;
1256 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
95808c2b 1257 tmpRecC.SetPtLeading(-1.);
d89b8229 1258 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1259 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1260 }// loop over random cones creation
1261
1262
1263 // loop over the reconstructed particles and add up the pT in the random cones
1264 // maybe better to loop over randomized particles not in the real jets...
1265 // but this by definition brings dow average energy in the whole event
1266 AliAODJet vTmpRanR(1,0,0,1);
1267 for(int i = 0; i < recParticles.GetEntries(); i++){
1268 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1269 if(fTCARandomConesOut){
1270 for(int ir = 0;ir < fNRandomCones;ir++){
1271 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1272 if(jC&&jC->DeltaR(vp)<fRparam){
1273 if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1274 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
95808c2b 1275 if(vp->Pt() > jC->GetPtLeading()) jC->SetPtLeading(vp->Pt());
d89b8229 1276 }
1277 }
1278 }// add up energy in cone
95808c2b 1279
d89b8229 1280 // the randomized input changes eta and phi, but keeps the p_T
1281 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1282 Double_t pTR = vp->Pt();
1283 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1284 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1285
1286 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1287 Double_t pZR = pTR/TMath::Tan(thetaR);
1288
1289 Double_t pXR = pTR * TMath::Cos(phiR);
1290 Double_t pYR = pTR * TMath::Sin(phiR);
1291 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1292 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1293 if(fTCARandomConesOutRan){
1294 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1295 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1296 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1297 if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1298 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
95808c2b 1299 if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt());
d89b8229 1300 }
1301 }
1302 }
1303 }
1304 }// loop over recparticles
1305
1306 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1307 if(fTCARandomConesOut){
1308 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
95808c2b 1309 // rescale the momentum vectors for the random cones
d89b8229 1310
1311 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1312 if(rC){
1313 Double_t etaC = rC->Eta();
1314 Double_t phiC = rC->Phi();
1315 // massless jet, unit vector
1316 pTC = rC->ChargedBgEnergy();
1317 if(pTC<=0)pTC = 0.001; // for almost empty events
1318 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1319 Double_t pZC = pTC/TMath::Tan(thetaC);
1320 Double_t pXC = pTC * TMath::Cos(phiC);
1321 Double_t pYC = pTC * TMath::Sin(phiC);
1322 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1323 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1324 rC->SetBgEnergy(0,0);
1325 rC->SetEffArea(jetArea,0);
1326 }
1327 }
1328 }
1329 if(fTCARandomConesOutRan){
1330 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1331 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1332 // same wit random
1333 if(rC){
1334 Double_t etaC = rC->Eta();
1335 Double_t phiC = rC->Phi();
1336 // massless jet, unit vector
1337 pTC = rC->ChargedBgEnergy();
1338 if(pTC<=0)pTC = 0.001;// for almost empty events
1339 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1340 Double_t pZC = pTC/TMath::Tan(thetaC);
1341 Double_t pXC = pTC * TMath::Cos(phiC);
1342 Double_t pYC = pTC * TMath::Sin(phiC);
1343 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1344 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1345 rC->SetBgEnergy(0,0);
1346 rC->SetEffArea(jetArea,0);
1347 }
1348 }
1349 }
1350 }// if(fNRandomCones
1351
1352 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1353
1354
1355
1356
1357
1358 if(fAODJetBackgroundOut){
1359 vector<fastjet::PseudoJet> jets2=sortedJets;
1360 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1361 Double_t bkg1=0;
1362 Double_t sigma1=0.;
1363 Double_t meanarea1=0.;
1364 Double_t bkg2=0;
1365 Double_t sigma2=0.;
1366 Double_t meanarea2=0.;
1367
1368 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1369 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1370
1371 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1372 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1373
1374 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1375 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1376 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1377 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1378
c0fdfe53 1379
1380 //Calculate rho with 4-vector area method for different quadrants with respect to the leading track in the event
1381 //exclude 2 leading jets from event
1382 //Quadrant 1: |phi_leadingTrack - phi_bkgCluster| < pi/2./2. - R (Near side to leading track)
1383 //Quadrant 2: pi/2 - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < pi/2 + (pi/2./2. - R) (Orthogonal to leading track)
1384 //Quadrant 3: pi - (pi/2/2 - R) < |phi_leadingTrack - phi_bkgCluster| < pi + (pi/2./2. - R) (Away side to leading track)
1385 //Quadrant 4: 3/2*pi - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < 3/2*pi + (pi/2./2. - R) (Orthogonal to leading track)
1386
1387 //Assuming R=0.2 for background clusters
1388
1389 Double_t bkg2Q[4] = {0.};
1390 Double_t sigma2Q[4] = {0.};
1391 Double_t meanarea2Q[4] = {0.};
1392
1393 //Get phi of leading track in event
1394 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1395 Float_t phiLeadingTrack = leading->Phi();
1396
1397 //Quadrant1 - near side
1398 phiMin = phiLeadingTrack - (TMath::Pi()/2./2. - 0.2);
1399 phiMax = phiLeadingTrack + (TMath::Pi()/2./2. - 0.2);
1400 fastjet::RangeDefinition rangeQ1(rapMin,rapMax, phiMin, phiMax);
1401 clustSeq.get_median_rho_and_sigma(jets2, rangeQ1, false, bkg2Q[0], sigma2Q[0], meanarea2Q[0], true);
1402
1403 //Quadrant2 - orthogonal
1404 phiMin = phiLeadingTrack + TMath::Pi()/2. - (TMath::Pi()/2./2. - 0.2);
1405 phiMax = phiLeadingTrack + TMath::Pi()/2. + (TMath::Pi()/2./2. - 0.2);
1406 fastjet::RangeDefinition rangeQ2(rapMin,rapMax, phiMin, phiMax);
1407 clustSeq.get_median_rho_and_sigma(jets2, rangeQ2, false, bkg2Q[1], sigma2Q[1], meanarea2Q[1], true);
1408
1409 //Quadrant3 - away side
1410 phiMin = phiLeadingTrack + TMath::Pi() - (TMath::Pi()/2./2. - 0.2);
1411 phiMax = phiLeadingTrack + TMath::Pi() + (TMath::Pi()/2./2. - 0.2);
1412 fastjet::RangeDefinition rangeQ3(rapMin,rapMax, phiMin, phiMax);
1413 clustSeq.get_median_rho_and_sigma(jets2, rangeQ3, false, bkg2Q[2], sigma2Q[2], meanarea2Q[2], true);
1414
1415 //Quadrant4 - othogonal
1416 phiMin = phiLeadingTrack + TMath::Pi()*3./2. - (TMath::Pi()/2./2. - 0.2);
1417 phiMax = phiLeadingTrack + TMath::Pi()*3./2. + (TMath::Pi()/2./2. - 0.2);
1418 fastjet::RangeDefinition rangeQ4(rapMin,rapMax, phiMin, phiMax);
1419 clustSeq.get_median_rho_and_sigma(jets2, rangeQ4, false, bkg2Q[3], sigma2Q[3], meanarea2Q[3], true);
1420
1421 fh3CentvsRhoLeadingTrackCorr->Fill(cent,bkg2Q[0],0.5);
1422 fh3CentvsRhoLeadingTrackCorr->Fill(cent,bkg2Q[1],1.5);
1423 fh3CentvsRhoLeadingTrackCorr->Fill(cent,bkg2Q[2],2.5);
1424 fh3CentvsRhoLeadingTrackCorr->Fill(cent,bkg2Q[3],3.5);
1425
1426 fh3CentvsSigmaLeadingTrackCorr->Fill(cent,sigma2Q[0],0.5);
1427 fh3CentvsSigmaLeadingTrackCorr->Fill(cent,sigma2Q[1],1.5);
1428 fh3CentvsSigmaLeadingTrackCorr->Fill(cent,sigma2Q[2],2.5);
1429 fh3CentvsSigmaLeadingTrackCorr->Fill(cent,sigma2Q[3],3.5);
1430
1431 fh3MultvsRhoLeadingTrackCorr->Fill(nCh,bkg2Q[0],0.5);
1432 fh3MultvsRhoLeadingTrackCorr->Fill(nCh,bkg2Q[1],1.5);
1433 fh3MultvsRhoLeadingTrackCorr->Fill(nCh,bkg2Q[2],2.5);
1434 fh3MultvsRhoLeadingTrackCorr->Fill(nCh,bkg2Q[3],3.5);
1435
1436 fh3MultvsSigmaLeadingTrackCorr->Fill(nCh,sigma2Q[0],0.5);
1437 fh3MultvsSigmaLeadingTrackCorr->Fill(nCh,sigma2Q[1],1.5);
1438 fh3MultvsSigmaLeadingTrackCorr->Fill(nCh,sigma2Q[2],2.5);
1439 fh3MultvsSigmaLeadingTrackCorr->Fill(nCh,sigma2Q[3],3.5);
1440
d89b8229 1441 }
1442 }
1443
1444
1445
1446
1447
1448 // fill track information
1449 Int_t nTrackOver = recParticles.GetSize();
1450 // do the same for tracks and jets
1451
1452 if(nTrackOver>0){
1453 TIterator *recIter = recParticles.MakeIterator();
1454 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1455 Float_t pt = tmpRec->Pt();
1456
1457 // Printf("Leading track p_t %3.3E",pt);
1458 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1459 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1460 while(pt<ptCut&&tmpRec){
1461 nTrackOver--;
1462 tmpRec = (AliVParticle*)(recIter->Next());
1463 if(tmpRec){
1464 pt = tmpRec->Pt();
1465 }
1466 }
1467 if(nTrackOver<=0)break;
1468 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1469 }
1470
1471 recIter->Reset();
1472 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1473 Float_t phi = leading->Phi();
1474 if(phi<0)phi+=TMath::Pi()*2.;
1475 Float_t eta = leading->Eta();
1476 pt = leading->Pt();
1477 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1478 Float_t tmpPt = tmpRec->Pt();
1479 Float_t tmpEta = tmpRec->Eta();
1480 fh1PtTracksRecIn->Fill(tmpPt);
1481 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1482 if(tmpRec==leading){
1483 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1484 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1485 continue;
1486 }
1487 // correlation
1488 Float_t tmpPhi = tmpRec->Phi();
1489
1490 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1491 Float_t dPhi = phi - tmpPhi;
1492 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1493 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1494 Float_t dEta = eta - tmpRec->Eta();
1495 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1496 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1497 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1498 }
1499 delete recIter;
1500 }
1501
1502 // find the random jets
1503
1504 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1505
1506 // fill the jet information from random track
1507 const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1508 const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1509
1510 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1511
1512 // loop over all jets an fill information, first on is the leading jet
1513
1514 Int_t nRecOverRan = inclusiveJetsRan.size();
1515 Int_t nRecRan = inclusiveJetsRan.size();
1516
1517 if(inclusiveJetsRan.size()>0){
1518 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1519 Float_t pt = leadingJet.Pt();
1520
1521 Int_t iCount = 0;
1522 TLorentzVector vecarearanb;
1523
1524 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1525 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1526 while(pt<ptCut&&iCount<nRecRan){
1527 nRecOverRan--;
1528 iCount++;
1529 if(iCount<nRecRan){
1530 pt = sortedJetsRan[iCount].perp();
1531 }
1532 }
1533 if(nRecOverRan<=0)break;
1534 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1535 }
1536 Float_t phi = leadingJet.Phi();
1537 if(phi<0)phi+=TMath::Pi()*2.;
1538 pt = leadingJet.Pt();
1539
1540 // correlation of leading jet with random tracks
1541
1542 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1543 {
1544 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1545 // correlation
1546 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1547 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1548 Float_t dPhi = phi - tmpPhi;
1549 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1550 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1551 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1552 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1553 }
1554
1555 Int_t nAodOutJetsRan = 0;
1556 AliAODJet *aodOutJetRan = 0;
1557 for(int j = 0; j < nRecRan;j++){
1558 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1559 Float_t tmpPt = tmpRec.Pt();
1560 fh1PtJetsRecInRan->Fill(tmpPt);
1561 // Fill Spectra with constituents
1562 const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1563 fh1NConstRecRan->Fill(constituents.size());
1564 fh2NConstPtRan->Fill(tmpPt,constituents.size());
1565 fh2PtNchRan->Fill(nCh,tmpPt);
1566 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1567
1568
1569 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1570 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1571 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1572 aodOutJetRan->GetRefTracks()->Clear();
1573 aodOutJetRan->SetEffArea(arearan,0);
1574 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1575 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1576 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1577
1578 }
1579
1580 // correlation
1581 Float_t tmpPhi = tmpRec.Phi();
1582 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1583
1584 if(j==0){
1585 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1586 fh1NConstLeadingRecRan->Fill(constituents.size());
1587 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1588 continue;
1589 }
1590 }
1591
1592
1593 if(fAODJetBackgroundOut){
1594 Double_t bkg3=0.;
1595 Double_t sigma3=0.;
1596 Double_t meanarea3=0.;
1597 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1598 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1599 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1600 /*
1601 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1602 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1603 */
1604 }
1605
1606
1607
1608 }
1609
1610
1611 // do the event selection if activated
1612 if(fJetTriggerPtCut>0){
1613 bool select = false;
1614 Float_t minPt = fJetTriggerPtCut;
1615 /*
1616 // hard coded for now ...
1617 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1618 if(cent<10)minPt = 50;
1619 else if(cent<30)minPt = 42;
1620 else if(cent<50)minPt = 28;
1621 else if(cent<80)minPt = 18;
1622 */
1623 float rho = 0;
1624 if(externalBackground)rho = externalBackground->GetBackground(2);
1625 if(fTCAJetsOut){
1626 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1627 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1628 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1629 if(ptSub>=minPt){
1630 select = true;
1631 break;
1632 }
1633 }
1634 }
1635
1636 if(select){
1637 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1638 fh1CentralitySelect->Fill(cent);
1639 fh1ZSelect->Fill(zVtx);
1640 aodH->SetFillAOD(kTRUE);
1641 }
1642 }
1643 if (fDebug > 2){
1644 if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1645 if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1646 if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1647 if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1648 }
1649 PostData(1, fHistList);
1650}
1651
1652void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1653{
f0659f11 1654 //
1655 // Terminate analysis
1656 //
006b2a30 1657 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1658
1659 if(fMomResH1Fit) delete fMomResH1Fit;
1660 if(fMomResH2Fit) delete fMomResH2Fit;
1661 if(fMomResH3Fit) delete fMomResH3Fit;
1662
d89b8229 1663}
1664
1665
1666Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1667
f0659f11 1668 //
1669 // get list of tracks/particles for different types
1670 //
1671
d89b8229 1672 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1673
1674 Int_t iCount = 0;
1675 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1676 if(type!=kTrackAODextraonly) {
1677 AliAODEvent *aod = 0;
1678 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1679 else aod = AODEvent();
1680 if(!aod){
1681 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1682 return iCount;
1683 }
d7396b04 1684
d89b8229 1685 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1686 AliAODTrack *tr = aod->GetTrack(it);
1687 Bool_t bGood = false;
1688 if(fFilterType == 0)bGood = true;
1689 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1690 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1691 if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1692 if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1693 continue;
1694 }
1695 if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1696 if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1697 continue;
1698 }
1699 if(tr->Pt()<fTrackPtCut){
1700 if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1701 continue;
1702 }
1703 if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1704 list->Add(tr);
1705 iCount++;
1706 }
1707 }
1708 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1709 AliAODEvent *aod = 0;
1710 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1711 else aod = AODEvent();
1712
1713 if(!aod){
1714 return iCount;
1715 }
1716 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1717 if(!aodExtraTracks)return iCount;
1718 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1719 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1720 if (!track) continue;
1721
1722 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1723 if(!trackAOD)continue;
1724 Bool_t bGood = false;
1725 if(fFilterType == 0)bGood = true;
1726 else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1727 else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1728 if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1729 if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1730 if(trackAOD->Pt()<fTrackPtCut) continue;
1731 list->Add(trackAOD);
1732 iCount++;
1733 }
1734 }
1735 }
1736 else if (type == kTrackKineAll||type == kTrackKineCharged){
1737 AliMCEvent* mcEvent = MCEvent();
1738 if(!mcEvent)return iCount;
1739 // we want to have alivpartilces so use get track
1740 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1741 if(!mcEvent->IsPhysicalPrimary(it))continue;
1742 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1743 if(type == kTrackKineAll){
1744 if(part->Pt()<fTrackPtCut)continue;
1745 list->Add(part);
1746 iCount++;
1747 }
1748 else if(type == kTrackKineCharged){
1749 if(part->Particle()->GetPDG()->Charge()==0)continue;
1750 if(part->Pt()<fTrackPtCut)continue;
1751 list->Add(part);
1752 iCount++;
1753 }
1754 }
1755 }
1756 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1757 AliAODEvent *aod = 0;
1758 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1759 else aod = AODEvent();
1760 if(!aod)return iCount;
1761 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1762 if(!tca)return iCount;
1763 for(int it = 0;it < tca->GetEntriesFast();++it){
1764 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1765 if(!part->IsPhysicalPrimary())continue;
1766 if(type == kTrackAODMCAll){
1767 if(part->Pt()<fTrackPtCut)continue;
1768 list->Add(part);
1769 iCount++;
1770 }
1771 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1772 if(part->Charge()==0)continue;
1773 if(part->Pt()<fTrackPtCut)continue;
1774 if(kTrackAODMCCharged){
1775 list->Add(part);
1776 }
1777 else {
1778 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1779 list->Add(part);
1780 }
1781 iCount++;
1782 }
1783 }
1784 }// AODMCparticle
1785 list->Sort();
1786 return iCount;
1787}
1788
fa7c34ba 1789void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() {
1790
575dee18 1791 TFile *f = TFile::Open(fPathTrPtResolution.Data());
fa7c34ba 1792
0e5bbbbc 1793 if(!f)return;
1794
fa7c34ba 1795 TProfile *fProfPtPtSigma1PtGlobSt = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
1796 TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
1797 TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
1798
1799 SetSmearResolution(kTRUE);
1800 SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD);
575dee18 1801
fa7c34ba 1802
1803}
1804
1805void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() {
1806
575dee18 1807 TFile *f = TFile::Open(fPathTrEfficiency.Data());
0e5bbbbc 1808 if(!f)return;
fa7c34ba 1809
1810 TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
1811 TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
1812 TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
1813
1814 SetDiceEfficiency(kTRUE);
1815
1816 if(fChangeEfficiencyFraction>0.) {
1817
1818 TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin");
1819
1820 for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) {
1821 Double_t content = hEffPosGlobSt->GetBinContent(bin);
1822 hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction);
1823 }
1824
1825 SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
1826
1827 }
1828 else
1829 SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
1830
fa7c34ba 1831}
1832
006b2a30 1833void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
1834
1835 //
1836 // set mom res profiles
1837 //
1838
575dee18 1839 if(fMomResH1) delete fMomResH1;
1840 if(fMomResH2) delete fMomResH2;
1841 if(fMomResH3) delete fMomResH3;
1842
1843 fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
1844 fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
1845 fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
1846
006b2a30 1847}
1848
1849void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
1850 //
1851 // set tracking efficiency histos
1852 //
1853
1854 fhEffH1 = (TH1*)h1->Clone("fhEffH1");
1855 fhEffH2 = (TH1*)h2->Clone("fhEffH2");
1856 fhEffH3 = (TH1*)h3->Clone("fhEffH3");
1857}
1858
1859Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
1860
1861 //
1862 // Get smearing on generated momentum
1863 //
1864
1865 //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
1866
1867 TProfile *fMomRes = 0x0;
1868 if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
1869 if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
1870 if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
1871
1872 if(!fMomRes) {
1873 return 0.;
1874 }
1875
1876
1877 Double_t smear = 0.;
1878
1879 if(pt>20.) {
1880 if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
1881 if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
1882 if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
1883 }
1884 else {
1885
1886 Int_t bin = fMomRes->FindBin(pt);
1887
1888 smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
1889
1890 }
1891
1892 if(fMomRes) delete fMomRes;
1893
1894 return smear;
1895}
1896
1897void AliAnalysisTaskJetCluster::FitMomentumResolution() {
1898 //
1899 // Fit linear function on momentum resolution at high pT
1900 //
1901
1902 if(!fMomResH1Fit && fMomResH1) {
1903 fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
1904 fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
1905 fMomResH1Fit ->SetRange(5.,100.);
1906 }
1907
1908 if(!fMomResH2Fit && fMomResH2) {
1909 fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
1910 fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
1911 fMomResH2Fit ->SetRange(5.,100.);
1912 }
1913
1914 if(!fMomResH3Fit && fMomResH3) {
1915 fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
1916 fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
1917 fMomResH3Fit ->SetRange(5.,100.);
1918 }
1919
1920}
1921
d89b8229 1922/*
1923Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1924 for(int i = 0; i < particles.GetEntries(); i++){
1925 AliVParticle *vp = (AliVParticle*)particles.At(i);
1926 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1927 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1928 jInp.set_user_index(i);
1929 inputParticles.push_back(jInp);
1930 }
1931
1932 return 0;
1933
1934}
1935*/