1 // **************************************
2 // Task used for the correction of determiantion of reconstructed jet spectra
3 // Compares input (gen) and output (rec) jets
4 // *******************************************
7 /**************************************************************************
8 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
10 * Author: The ALICE Off-line Project. *
11 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
26 #include <TInterpreter.h>
28 #include <TRefArray.h>
37 #include <TLorentzVector.h>
38 #include <TClonesArray.h>
39 #include "TDatabasePDG.h"
42 #include "AliAnalysisTaskJetCluster.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 "AliAODMCParticle.h"
54 #include "AliMCEventHandler.h"
55 #include "AliMCEvent.h"
57 #include "AliGenPythiaEventHeader.h"
58 #include "AliJetKineReaderHeader.h"
59 #include "AliGenCocktailEventHeader.h"
60 #include "AliInputEventHandler.h"
61 #include "AliAODJetEventBackground.h"
63 #include "fastjet/PseudoJet.hh"
64 #include "fastjet/ClusterSequenceArea.hh"
65 #include "fastjet/AreaDefinition.hh"
66 #include "fastjet/JetDefinition.hh"
67 // get info on how fastjet was configured
68 #include "fastjet/config.h"
72 ClassImp(AliAnalysisTaskJetCluster)
74 AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
82 if(fTCAJetsOut)fTCAJetsOut->Delete();
85 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
86 delete fTCAJetsOutRan;
88 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
89 delete fTCARandomConesOut;
91 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
92 delete fTCARandomConesOutRan;
94 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
95 delete fAODJetBackgroundOut;
98 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster():
103 fUseAODTrackInput(kFALSE),
104 fUseAODMCInput(kFALSE),
105 fUseBackgroundCalc(kFALSE),
106 fEventSelection(kFALSE),
107 fRequireVZEROAC(kFALSE),
108 fRequireTZEROvtx(kFALSE),
110 fFilterMaskBestPt(0),
113 fTrackTypeRec(kTrackUndef),
114 fTrackTypeGen(kTrackUndef),
116 fNSkipLeadingCone(0),
120 fTrackEtaWindow(0.9),
123 fJetOutputMinPt(0.150),
124 fMaxTrackPtInJet(100.),
130 fStoreRhoLeadingTrackCorr(kFALSE),
132 fBackgroundBranch(""),
143 fUseTrPtResolutionSmearing(kFALSE),
144 fUseDiceEfficiency(kFALSE),
145 fDiceEfficiencyMinPt(-1.),
146 fUseTrPtResolutionFromOADB(kFALSE),
147 fUseTrEfficiencyFromOADB(kFALSE),
148 fPathTrPtResolution(""),
149 fPathTrEfficiency(""),
150 fChangeEfficiencyFraction(0.),
151 fEfficiencyFixed(1.),
153 fAlgorithm(fastjet::kt_algorithm),
154 fStrategy(fastjet::Best),
155 fRecombScheme(fastjet::BIpt_scheme),
156 fAreaType(fastjet::active_area),
158 fActiveAreaRepeats(1),
162 fTCARandomConesOut(0x0),
163 fTCARandomConesOutRan(0x0),
164 fAODJetBackgroundOut(0x0),
170 fh1PtHardTrials(0x0),
173 fh1NConstLeadingRec(0x0),
175 fh1PtJetsLeadingRecIn(0x0),
176 fh1PtJetConstRec(0x0),
177 fh1PtJetConstLeadingRec(0x0),
178 fh1PtTracksRecIn(0x0),
179 fh1PtTracksLeadingRecIn(0x0),
181 fh1NConstRecRan(0x0),
182 fh1PtJetsLeadingRecInRan(0x0),
183 fh1NConstLeadingRecRan(0x0),
184 fh1PtJetsRecInRan(0x0),
185 fh1PtTracksGenIn(0x0),
187 fh1CentralityPhySel(0x0),
189 fh1CentralitySelect(0x0),
194 fh2NRecTracksPt(0x0),
196 fh2NConstLeadingPt(0x0),
198 fh2LeadingJetPhiEta(0x0),
200 fh2LeadingJetEtaPt(0x0),
202 fh2LeadingTrackEtaPt(0x0),
203 fh2JetsLeadingPhiEta(0x0),
204 fh2JetsLeadingPhiPt(0x0),
205 fh2TracksLeadingPhiEta(0x0),
206 fh2TracksLeadingPhiPt(0x0),
207 fh2TracksLeadingJetPhiPt(0x0),
208 fh2JetsLeadingPhiPtW(0x0),
209 fh2TracksLeadingPhiPtW(0x0),
210 fh2TracksLeadingJetPhiPtW(0x0),
211 fh2NRecJetsPtRan(0x0),
213 fh2NConstLeadingPtRan(0x0),
218 fh2TracksLeadingJetPhiPtRan(0x0),
219 fh2TracksLeadingJetPhiPtWRan(0x0),
220 fh3CentvsRhoLeadingTrackPt(0x0),
221 fh3CentvsSigmaLeadingTrackPt(0x0),
222 fh3MultvsRhoLeadingTrackPt(0x0),
223 fh3MultvsSigmaLeadingTrackPt(0x0),
224 fh3CentvsRhoLeadingTrackPtQ1(0x0),
225 fh3CentvsRhoLeadingTrackPtQ2(0x0),
226 fh3CentvsRhoLeadingTrackPtQ3(0x0),
227 fh3CentvsRhoLeadingTrackPtQ4(0x0),
228 fh3CentvsSigmaLeadingTrackPtQ1(0x0),
229 fh3CentvsSigmaLeadingTrackPtQ2(0x0),
230 fh3CentvsSigmaLeadingTrackPtQ3(0x0),
231 fh3CentvsSigmaLeadingTrackPtQ4(0x0),
232 fh3MultvsRhoLeadingTrackPtQ1(0x0),
233 fh3MultvsRhoLeadingTrackPtQ2(0x0),
234 fh3MultvsRhoLeadingTrackPtQ3(0x0),
235 fh3MultvsRhoLeadingTrackPtQ4(0x0),
236 fh3MultvsSigmaLeadingTrackPtQ1(0x0),
237 fh3MultvsSigmaLeadingTrackPtQ2(0x0),
238 fh3MultvsSigmaLeadingTrackPtQ3(0x0),
239 fh3MultvsSigmaLeadingTrackPtQ4(0x0),
240 fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0),
241 fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0),
242 fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0),
243 fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0),
244 fh2PtGenPtSmeared(0x0),
246 fp1PtResolution(0x0),
253 for(int i = 0;i<3;i++){
254 fh1BiARandomCones[i] = 0;
255 fh1BiARandomConesRan[i] = 0;
257 for(int i = 0;i<kMaxCent;i++){
258 fh2JetsLeadingPhiPtC[i] = 0;
259 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
260 fh2TracksLeadingJetPhiPtC[i] = 0;
261 fh2TracksLeadingJetPhiPtWC[i] = 0;
265 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
266 AliAnalysisTaskSE(name),
270 fUseAODTrackInput(kFALSE),
271 fUseAODMCInput(kFALSE),
272 fUseBackgroundCalc(kFALSE),
273 fEventSelection(kFALSE),
274 fRequireVZEROAC(kFALSE),
275 fRequireTZEROvtx(kFALSE),
277 fFilterMaskBestPt(0),
280 fTrackTypeRec(kTrackUndef),
281 fTrackTypeGen(kTrackUndef),
283 fNSkipLeadingCone(0),
287 fTrackEtaWindow(0.9),
290 fJetOutputMinPt(0.150),
291 fMaxTrackPtInJet(100.),
297 fStoreRhoLeadingTrackCorr(kFALSE),
299 fBackgroundBranch(""),
310 fUseTrPtResolutionSmearing(kFALSE),
311 fUseDiceEfficiency(kFALSE),
312 fDiceEfficiencyMinPt(-1.),
313 fUseTrPtResolutionFromOADB(kFALSE),
314 fUseTrEfficiencyFromOADB(kFALSE),
315 fPathTrPtResolution(""),
316 fPathTrEfficiency(""),
317 fChangeEfficiencyFraction(0.),
318 fEfficiencyFixed(1.),
320 fAlgorithm(fastjet::kt_algorithm),
321 fStrategy(fastjet::Best),
322 fRecombScheme(fastjet::BIpt_scheme),
323 fAreaType(fastjet::active_area),
325 fActiveAreaRepeats(1),
329 fTCARandomConesOut(0x0),
330 fTCARandomConesOutRan(0x0),
331 fAODJetBackgroundOut(0x0),
337 fh1PtHardTrials(0x0),
340 fh1NConstLeadingRec(0x0),
342 fh1PtJetsLeadingRecIn(0x0),
343 fh1PtJetConstRec(0x0),
344 fh1PtJetConstLeadingRec(0x0),
345 fh1PtTracksRecIn(0x0),
346 fh1PtTracksLeadingRecIn(0x0),
348 fh1NConstRecRan(0x0),
349 fh1PtJetsLeadingRecInRan(0x0),
350 fh1NConstLeadingRecRan(0x0),
351 fh1PtJetsRecInRan(0x0),
352 fh1PtTracksGenIn(0x0),
354 fh1CentralityPhySel(0x0),
356 fh1CentralitySelect(0x0),
361 fh2NRecTracksPt(0x0),
363 fh2NConstLeadingPt(0x0),
365 fh2LeadingJetPhiEta(0x0),
367 fh2LeadingJetEtaPt(0x0),
369 fh2LeadingTrackEtaPt(0x0),
370 fh2JetsLeadingPhiEta(0x0),
371 fh2JetsLeadingPhiPt(0x0),
372 fh2TracksLeadingPhiEta(0x0),
373 fh2TracksLeadingPhiPt(0x0),
374 fh2TracksLeadingJetPhiPt(0x0),
375 fh2JetsLeadingPhiPtW(0x0),
376 fh2TracksLeadingPhiPtW(0x0),
377 fh2TracksLeadingJetPhiPtW(0x0),
378 fh2NRecJetsPtRan(0x0),
380 fh2NConstLeadingPtRan(0x0),
385 fh2TracksLeadingJetPhiPtRan(0x0),
386 fh2TracksLeadingJetPhiPtWRan(0x0),
387 fh3CentvsRhoLeadingTrackPt(0x0),
388 fh3CentvsSigmaLeadingTrackPt(0x0),
389 fh3MultvsRhoLeadingTrackPt(0x0),
390 fh3MultvsSigmaLeadingTrackPt(0x0),
391 fh3CentvsRhoLeadingTrackPtQ1(0x0),
392 fh3CentvsRhoLeadingTrackPtQ2(0x0),
393 fh3CentvsRhoLeadingTrackPtQ3(0x0),
394 fh3CentvsRhoLeadingTrackPtQ4(0x0),
395 fh3CentvsSigmaLeadingTrackPtQ1(0x0),
396 fh3CentvsSigmaLeadingTrackPtQ2(0x0),
397 fh3CentvsSigmaLeadingTrackPtQ3(0x0),
398 fh3CentvsSigmaLeadingTrackPtQ4(0x0),
399 fh3MultvsRhoLeadingTrackPtQ1(0x0),
400 fh3MultvsRhoLeadingTrackPtQ2(0x0),
401 fh3MultvsRhoLeadingTrackPtQ3(0x0),
402 fh3MultvsRhoLeadingTrackPtQ4(0x0),
403 fh3MultvsSigmaLeadingTrackPtQ1(0x0),
404 fh3MultvsSigmaLeadingTrackPtQ2(0x0),
405 fh3MultvsSigmaLeadingTrackPtQ3(0x0),
406 fh3MultvsSigmaLeadingTrackPtQ4(0x0),
407 fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0),
408 fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0),
409 fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0),
410 fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0),
411 fh2PtGenPtSmeared(0x0),
413 fp1PtResolution(0x0),
420 for(int i = 0;i<3;i++){
421 fh1BiARandomCones[i] = 0;
422 fh1BiARandomConesRan[i] = 0;
424 for(int i = 0;i<kMaxCent;i++){
425 fh2JetsLeadingPhiPtC[i] = 0;
426 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
427 fh2TracksLeadingJetPhiPtC[i] = 0;
428 fh2TracksLeadingJetPhiPtWC[i] = 0;
430 DefineOutput(1, TList::Class());
435 Bool_t AliAnalysisTaskJetCluster::Notify()
438 // Implemented Notify() to read the cross sections
439 // and number of trials from pyxsec.root
444 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
448 // Create the output container
451 fRandom = new TRandom3(0);
457 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
461 if(fNonStdBranch.Length()!=0)
463 // only create the output branch if we have a name
464 // Create a new branch for jets...
465 // -> cleared in the UserExec....
466 // here we can also have the case that the brnaches are written to a separate file
469 fTCAJetsOut = new TClonesArray("AliAODJet", 0);
470 fTCAJetsOut->SetName(fNonStdBranch.Data());
471 AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
474 if(fJetTypes&kJetRan){
475 fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
476 fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
477 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
478 if( fEfficiencyFixed < 1.)
479 fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dEffFixed%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.)));
481 fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dFr%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
483 AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
486 if(fUseBackgroundCalc){
487 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
488 fAODJetBackgroundOut = new AliAODJetEventBackground();
489 fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
490 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
491 if( fEfficiencyFixed < 1.)
492 fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dEffFixed%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.)));
494 fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dFr%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
496 AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
499 // create the branch for the random cones with the same R
500 TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
501 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
502 if( fEfficiencyFixed < 1.)
503 cName = Form("%sDetector%d%dEffFixed%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.),fNSkipLeadingCone);
505 cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone);
509 if(!AODEvent()->FindListObject(cName.Data())){
510 fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
511 fTCARandomConesOut->SetName(cName.Data());
512 AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
515 // create the branch with the random for the random cones on the random event
516 if(fJetTypes&kRCRan){
517 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
518 if(!AODEvent()->FindListObject(cName.Data())){
519 fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
520 fTCARandomConesOutRan->SetName(cName.Data());
521 AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
526 if(fNonStdFile.Length()!=0){
528 // case that we have an AOD extension we need to fetch the jets from the extended output
529 // we identify the extension aod event by looking for the branchname
530 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
531 // case that we have an AOD extension we need can fetch the background maybe from the extended output
532 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
537 if(!fHistList)fHistList = new TList();
538 fHistList->SetOwner();
539 PostData(1, fHistList); // post data in any case once
541 Bool_t oldStatus = TH1::AddDirectoryStatus();
542 TH1::AddDirectory(kFALSE);
547 const Int_t nBinPt = 100;
548 Double_t binLimitsPt[nBinPt+1];
549 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
551 binLimitsPt[iPt] = 0.0;
554 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
558 const Int_t nBinPhi = 90;
559 Double_t binLimitsPhi[nBinPhi+1];
560 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
562 binLimitsPhi[iPhi] = -1.*TMath::Pi();
565 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
571 const Int_t nBinEta = 40;
572 Double_t binLimitsEta[nBinEta+1];
573 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
575 binLimitsEta[iEta] = -2.0;
578 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
582 const int nChMax = 5000;
584 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
585 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
587 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
588 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
591 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
592 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
594 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
595 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
596 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
597 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
600 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
601 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
602 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
604 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
605 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
606 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
607 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
608 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
609 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
610 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
611 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
612 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
613 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
615 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
616 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
617 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
619 fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
620 fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
621 fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
623 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
624 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
625 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
629 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
630 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
631 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
632 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
634 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
635 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);
636 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);
637 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);
641 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
642 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
643 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
644 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
646 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
647 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
648 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
649 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
651 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
652 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
653 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
654 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
658 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
659 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
660 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
661 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
662 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
663 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
664 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
665 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
666 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
667 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
668 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
669 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
671 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
672 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
673 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
674 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
676 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
677 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
678 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
679 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
681 if(fStoreRhoLeadingTrackCorr) {
682 fh3CentvsRhoLeadingTrackPt = new TH3F("fh3CentvsRhoLeadingTrackPt","centrality vs background density full event; centrality; #rho", 50,0.,100.,500,0.,250.,100,0.,100.);
683 fh3CentvsSigmaLeadingTrackPt = new TH3F("fh3CentvsSigmaLeadingTrackPt","centrality vs sigma full event; centrality; #sigma(#rho)", 50,0.,100.,50,0.,50.,100,0.,100.);
684 fh3MultvsRhoLeadingTrackPt = new TH3F("fh3MultvsRhoLeadingTrackPt","multiplicity vs background density full event; multiplicity; #rho", 100,0.,5000.,500,0.,250.,100,0.,100.);
685 fh3MultvsSigmaLeadingTrackPt = new TH3F("fh3MultvsSigmaLeadingTrackPt","multiplicity vs sigma full event; multiplicity; #sigma(#rho)", 100,0.,5000.,50,0.,50.,100,0.,100.);
688 fh3CentvsRhoLeadingTrackPtQ1 = new TH3F("fh3CentvsRhoLeadingTrackPtQ1","centrality vs background density Q1; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
689 fh3CentvsRhoLeadingTrackPtQ2 = new TH3F("fh3CentvsRhoLeadingTrackPtQ2","centrality vs background density Q2; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
690 fh3CentvsRhoLeadingTrackPtQ3 = new TH3F("fh3CentvsRhoLeadingTrackPtQ3","centrality vs background density Q3; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
691 fh3CentvsRhoLeadingTrackPtQ4 = new TH3F("fh3CentvsRhoLeadingTrackPtQ4","centrality vs background density Q4; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
693 fh3CentvsSigmaLeadingTrackPtQ1 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ1","centrality vs background density Q1; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
694 fh3CentvsSigmaLeadingTrackPtQ2 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ2","centrality vs background density Q2; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
695 fh3CentvsSigmaLeadingTrackPtQ3 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ3","centrality vs background density Q3; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
696 fh3CentvsSigmaLeadingTrackPtQ4 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ4","centrality vs background density Q4; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
698 fh3MultvsRhoLeadingTrackPtQ1 = new TH3F("fh3MultvsRhoLeadingTrackPtQ1","multiplicity vs background density Q1; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
699 fh3MultvsRhoLeadingTrackPtQ2 = new TH3F("fh3MultvsRhoLeadingTrackPtQ2","multiplicity vs background density Q2; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
700 fh3MultvsRhoLeadingTrackPtQ3 = new TH3F("fh3MultvsRhoLeadingTrackPtQ3","multiplicity vs background density Q3; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
701 fh3MultvsRhoLeadingTrackPtQ4 = new TH3F("fh3MultvsRhoLeadingTrackPtQ4","multiplicity vs background density Q4; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
703 fh3MultvsSigmaLeadingTrackPtQ1 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ1","multiplicity vs background density Q1; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
704 fh3MultvsSigmaLeadingTrackPtQ2 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ2","multiplicity vs background density Q2; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
705 fh3MultvsSigmaLeadingTrackPtQ3 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ3","multiplicity vs background density Q3; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
706 fh3MultvsSigmaLeadingTrackPtQ4 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ4","multiplicity vs background density Q4; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
709 fh3CentvsDeltaRhoLeadingTrackPtQ1 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ1","centrality vs background density Q1; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.);
710 fh3CentvsDeltaRhoLeadingTrackPtQ2 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ2","centrality vs background density Q2; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.);
711 fh3CentvsDeltaRhoLeadingTrackPtQ3 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ3","centrality vs background density Q3; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.);
712 fh3CentvsDeltaRhoLeadingTrackPtQ4 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ4","centrality vs background density Q4; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.);
714 fHistList->Add(fh3CentvsRhoLeadingTrackPt);
715 fHistList->Add(fh3CentvsSigmaLeadingTrackPt);
716 fHistList->Add(fh3MultvsRhoLeadingTrackPt);
717 fHistList->Add(fh3MultvsSigmaLeadingTrackPt);
719 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ1);
720 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ2);
721 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ3);
722 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ4);
724 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ1);
725 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ2);
726 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ3);
727 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ4);
729 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ1);
730 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ2);
731 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ3);
732 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ4);
734 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ1);
735 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ2);
736 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ3);
737 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ4);
739 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ1);
740 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ2);
741 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ3);
742 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ4);
746 //Detector level effects histos
747 fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
749 fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
750 fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
752 fHistList->Add(fh2PtGenPtSmeared);
753 fHistList->Add(fp1Efficiency);
754 fHistList->Add(fp1PtResolution);
756 if(fNRandomCones>0&&fUseBackgroundCalc){
757 for(int i = 0;i<3;i++){
758 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
759 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
763 for(int i = 0;i < kMaxCent;i++){
764 fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
765 fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
766 fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
767 fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
770 const Int_t saveLevel = 3; // large save level more histos
772 fHistList->Add(fh1Xsec);
773 fHistList->Add(fh1Trials);
775 fHistList->Add(fh1NJetsRec);
776 fHistList->Add(fh1NConstRec);
777 fHistList->Add(fh1NConstLeadingRec);
778 fHistList->Add(fh1PtJetsRecIn);
779 fHistList->Add(fh1PtJetsLeadingRecIn);
780 fHistList->Add(fh1PtTracksRecIn);
781 fHistList->Add(fh1PtTracksLeadingRecIn);
782 fHistList->Add(fh1PtJetConstRec);
783 fHistList->Add(fh1PtJetConstLeadingRec);
784 fHistList->Add(fh1NJetsRecRan);
785 fHistList->Add(fh1NConstRecRan);
786 fHistList->Add(fh1PtJetsLeadingRecInRan);
787 fHistList->Add(fh1NConstLeadingRecRan);
788 fHistList->Add(fh1PtJetsRecInRan);
789 fHistList->Add(fh1Nch);
790 fHistList->Add(fh1Centrality);
791 fHistList->Add(fh1CentralitySelect);
792 fHistList->Add(fh1CentralityPhySel);
793 fHistList->Add(fh1Z);
794 fHistList->Add(fh1ZSelect);
795 fHistList->Add(fh1ZPhySel);
796 if(fNRandomCones>0&&fUseBackgroundCalc){
797 for(int i = 0;i<3;i++){
798 fHistList->Add(fh1BiARandomCones[i]);
799 fHistList->Add(fh1BiARandomConesRan[i]);
802 for(int i = 0;i < kMaxCent;i++){
803 fHistList->Add(fh2JetsLeadingPhiPtC[i]);
804 fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
805 fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
806 fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
809 fHistList->Add(fh2NRecJetsPt);
810 fHistList->Add(fh2NRecTracksPt);
811 fHistList->Add(fh2NConstPt);
812 fHistList->Add(fh2NConstLeadingPt);
813 fHistList->Add(fh2PtNch);
814 fHistList->Add(fh2PtNchRan);
815 fHistList->Add(fh2PtNchN);
816 fHistList->Add(fh2PtNchNRan);
817 fHistList->Add(fh2JetPhiEta);
818 fHistList->Add(fh2LeadingJetPhiEta);
819 fHistList->Add(fh2JetEtaPt);
820 fHistList->Add(fh2LeadingJetEtaPt);
821 fHistList->Add(fh2TrackEtaPt);
822 fHistList->Add(fh2LeadingTrackEtaPt);
823 fHistList->Add(fh2JetsLeadingPhiEta );
824 fHistList->Add(fh2JetsLeadingPhiPt);
825 fHistList->Add(fh2TracksLeadingPhiEta);
826 fHistList->Add(fh2TracksLeadingPhiPt);
827 fHistList->Add(fh2TracksLeadingJetPhiPt);
828 fHistList->Add(fh2JetsLeadingPhiPtW);
829 fHistList->Add(fh2TracksLeadingPhiPtW);
830 fHistList->Add(fh2TracksLeadingJetPhiPtW);
831 fHistList->Add(fh2NRecJetsPtRan);
832 fHistList->Add(fh2NConstPtRan);
833 fHistList->Add(fh2NConstLeadingPtRan);
834 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
835 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
838 // =========== Switch on Sumw2 for all histos ===========
839 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
840 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
845 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
848 TH1::AddDirectory(oldStatus);
851 void AliAnalysisTaskJetCluster::LocalInit()
857 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
859 if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB();
860 if(fUseTrEfficiencyFromOADB) LoadTrEfficiencyRootFileFromOADB();
863 FitMomentumResolution();
867 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
870 // handle and reset the output jet branch
872 if(fTCAJetsOut)fTCAJetsOut->Delete();
873 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
874 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
875 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
876 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
878 AliAODJetEventBackground* externalBackground = 0;
879 if(!externalBackground&&fBackgroundBranch.Length()){
880 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
881 if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
882 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
885 // Execute analysis for current event
887 AliESDEvent *fESD = 0;
888 if(fUseAODTrackInput){
889 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
891 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
897 // assume that the AOD is in the general output...
900 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
904 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
908 //Check if information is provided detector level effects
909 if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE;
910 if(!fhEffH1 || !fhEffH2 || !fhEffH3 ) fUseDiceEfficiency = kFALSE;
911 if( fEfficiencyFixed < 1. ) fUseDiceEfficiency = kTRUE;
913 Bool_t selectEvent = false;
914 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
920 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
921 TString vtxTitle(vtxAOD->GetTitle());
922 zVtx = vtxAOD->GetZ();
924 cent = fAOD->GetHeader()->GetCentrality();
925 if(cent<10)cenClass = 0;
926 else if(cent<30)cenClass = 1;
927 else if(cent<50)cenClass = 2;
928 else if(cent<80)cenClass = 3;
929 if(physicsSelection){
930 fh1CentralityPhySel->Fill(cent);
931 fh1ZPhySel->Fill(zVtx);
935 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
936 Float_t yvtx = vtxAOD->GetY();
937 Float_t xvtx = vtxAOD->GetX();
938 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
939 if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
940 if(physicsSelection){
946 if(cent<fCentCutLo||cent>fCentCutUp){
958 const AliAODVZERO *vzero = fAOD->GetVZEROData();
960 if((vzero->GetTriggerChargeA()>0)&&(vzero->GetTriggerChargeC()>0)){
965 const AliAODTZERO *tzero = fAOD->GetTZEROData();
967 if(TMath::Abs(tzero->GetT0VertexRaw())<100){
972 if(fRequireVZEROAC&&fRequireTZEROvtx)selectEvent = selectEvent&&V0&&T0;
973 else if(fRequireTZEROvtx)selectEvent = selectEvent&&T0;
974 else if(fRequireVZEROAC)selectEvent = selectEvent&&V0;
978 PostData(1, fHistList);
981 fh1Centrality->Fill(cent);
983 fh1Trials->Fill("#sum{ntrials}",1);
986 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
988 // ==== General variables needed
992 // we simply fetch the tracks/mc particles as a list of AliVParticles
997 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
998 Float_t nCh = recParticles.GetEntries();
1000 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
1001 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
1002 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
1004 // find the jets....
1006 vector<fastjet::PseudoJet> inputParticlesRec;
1007 vector<fastjet::PseudoJet> inputParticlesRecRan;
1009 // Generate the random cones
1011 AliAODJet vTmpRan(1,0,0,1);
1012 for(int i = 0; i < recParticles.GetEntries(); i++){
1013 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1015 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1016 // we take total momentum here
1018 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
1019 //Add particles to fastjet in case we are not running toy model
1020 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
1021 jInp.set_user_index(i);
1022 inputParticlesRec.push_back(jInp);
1024 else if(fUseDiceEfficiency) {
1026 // Dice to decide if particle is kept or not - toy model for efficiency
1028 Double_t sumEff = 0.;
1030 Double_t eff[3] = {0.};
1033 Double_t rnd = fRandom->Uniform(1.);
1034 if( fEfficiencyFixed<1. ) {
1035 sumEff = fEfficiencyFixed;
1039 Double_t pTtmp = pT;
1040 if(pT>10.) pTtmp = 10.;
1041 Double_t eff1 = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
1042 Double_t eff2 = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
1043 Double_t eff3 = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));
1045 //Sort efficiencies from large to small
1046 if(eff1>eff2 && eff1>eff3) {
1061 else if(eff2>eff1 && eff2>eff3) {
1076 else if(eff3>eff1 && eff3>eff2) {
1092 sumEff = eff[0]+eff[1]+eff[2];
1094 fp1Efficiency->Fill(vp->Pt(),sumEff);
1095 if(rnd>sumEff && pT > fDiceEfficiencyMinPt) continue;
1097 if(fUseTrPtResolutionSmearing) {
1098 //Smear momentum of generated particle
1099 Double_t smear = 1.;
1100 //Select hybrid track category
1102 smear = GetMomentumSmearing(cat[2],pT);
1103 else if(rnd<=(eff[2]+eff[1]))
1104 smear = GetMomentumSmearing(cat[1],pT);
1106 smear = GetMomentumSmearing(cat[0],pT);
1108 fp1PtResolution->Fill(vp->Pt(),smear);
1110 Double_t sigma = vp->Pt()*smear;
1111 Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
1113 Double_t phi = vp->Phi();
1114 Double_t theta = 2.*TMath::ATan(TMath::Exp(-1.*(vp->Eta())));
1115 Double_t pX = pTrec * TMath::Cos(phi);
1116 Double_t pY = pTrec * TMath::Sin(phi);
1117 Double_t pZ = pTrec/TMath::Tan(theta);
1118 Double_t p=TMath::Sqrt(pTrec*pTrec+pZ*pZ);
1120 fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);
1122 fastjet::PseudoJet jInp(pX,pY,pZ,p);
1123 jInp.set_user_index(i);
1124 inputParticlesRec.push_back(jInp);
1128 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
1129 jInp.set_user_index(i);
1130 inputParticlesRec.push_back(jInp);
1136 // the randomized input changes eta and phi, but keeps the p_T
1137 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1138 Double_t pT = vp->Pt();
1139 Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
1140 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
1142 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
1143 Double_t pZ = pT/TMath::Tan(theta);
1145 Double_t pX = pT * TMath::Cos(phi);
1146 Double_t pY = pT * TMath::Sin(phi);
1147 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
1148 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
1150 jInpRan.set_user_index(i);
1151 inputParticlesRecRan.push_back(jInpRan);
1152 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
1155 // fill the tref array, only needed when we write out jets
1158 fRef->Delete(); // make sure to delete before placement new...
1159 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
1160 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
1163 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ...
1167 if(inputParticlesRec.size()==0){
1168 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
1169 PostData(1, fHistList);
1174 // employ setters for these...
1177 // now create the object that holds info about ghosts
1179 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
1180 // reduce CPU time...
1182 fActiveAreaRepeats = 0;
1186 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
1187 fastjet::AreaType areaType = fastjet::active_area;
1188 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
1189 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
1190 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
1192 //range where to compute background
1193 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
1195 phiMax = 2*TMath::Pi();
1196 rapMax = fGhostEtamax - fRparam;
1197 rapMin = - fGhostEtamax + fRparam;
1198 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
1201 const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
1202 const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
1205 fh1NJetsRec->Fill(sortedJets.size());
1207 // loop over all jets an fill information, first on is the leading jet
1209 Int_t nRecOver = inclusiveJets.size();
1210 Int_t nRec = inclusiveJets.size();
1211 if(inclusiveJets.size()>0){
1212 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
1213 Double_t area = clustSeq.area(sortedJets[0]);
1214 leadingJet.SetEffArea(area,0);
1215 Float_t pt = leadingJet.Pt();
1216 Int_t nAodOutJets = 0;
1217 Int_t nAodOutTracks = 0;
1218 AliAODJet *aodOutJet = 0;
1221 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
1222 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
1223 while(pt<ptCut&&iCount<nRec){
1227 pt = sortedJets[iCount].perp();
1230 if(nRecOver<=0)break;
1231 fh2NRecJetsPt->Fill(ptCut,nRecOver);
1233 Float_t phi = leadingJet.Phi();
1234 if(phi<0)phi+=TMath::Pi()*2.;
1235 Float_t eta = leadingJet.Eta();
1237 if(externalBackground){
1238 // carefull has to be filled in a task before
1239 // todo, ReArrange to the botom
1240 pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
1242 pt = leadingJet.Pt() - pTback;
1243 // correlation of leading jet with tracks
1244 TIterator *recIter = recParticles.MakeIterator();
1246 AliVParticle *tmpRecTrack = 0;
1247 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
1248 Float_t tmpPt = tmpRecTrack->Pt();
1250 Float_t tmpPhi = tmpRecTrack->Phi();
1251 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1252 Float_t dPhi = phi - tmpPhi;
1253 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1254 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1255 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
1256 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
1258 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
1259 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1265 TLorentzVector vecareab;
1266 for(int j = 0; j < nRec;j++){
1267 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
1270 Float_t tmpPt = tmpRec.Pt();
1272 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
1273 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
1274 aodOutJet->GetRefTracks()->Clear();
1275 Double_t area1 = clustSeq.area(sortedJets[j]);
1276 aodOutJet->SetEffArea(area1,0);
1277 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
1278 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
1279 aodOutJet->SetVectorAreaCharged(&vecareab);
1283 Float_t tmpPtBack = 0;
1284 if(externalBackground){
1285 // carefull has to be filled in a task before
1286 // todo, ReArrange to the botom
1287 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
1289 tmpPt = tmpPt - tmpPtBack;
1290 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
1292 fh1PtJetsRecIn->Fill(tmpPt);
1293 // Fill Spectra with constituentsemacs
1294 const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
1296 fh1NConstRec->Fill(constituents.size());
1297 fh2PtNch->Fill(nCh,tmpPt);
1298 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
1299 fh2NConstPt->Fill(tmpPt,constituents.size());
1300 // loop over constiutents and fill spectrum
1302 AliVParticle *partLead = 0;
1303 Float_t ptLead = -1;
1305 for(unsigned int ic = 0; ic < constituents.size();ic++){
1306 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1308 fh1PtJetConstRec->Fill(part->Pt());
1310 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1311 if(part->Pt()>fMaxTrackPtInJet){
1312 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1315 if(part->Pt()>ptLead){
1316 ptLead = part->Pt();
1319 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1322 AliAODTrack *aodT = 0;
1325 //set pT of leading constituent of jet
1326 aodOutJet->SetPtLeading(partLead->Pt());
1327 aodT = dynamic_cast<AliAODTrack*>(partLead);
1329 if(aodT->TestFilterBit(fFilterMaskBestPt)){
1330 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
1337 Float_t tmpPhi = tmpRec.Phi();
1338 Float_t tmpEta = tmpRec.Eta();
1339 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1341 fh1PtJetsLeadingRecIn->Fill(tmpPt);
1342 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1343 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1344 fh1NConstLeadingRec->Fill(constituents.size());
1345 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1348 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1349 fh2JetEtaPt->Fill(tmpEta,tmpPt);
1350 Float_t dPhi = phi - tmpPhi;
1351 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1352 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1353 Float_t dEta = eta - tmpRec.Eta();
1354 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1355 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1357 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1358 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1360 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1361 }// loop over reconstructed jets
1366 // Add the random cones...
1367 if(fNRandomCones>0&&fTCARandomConesOut){
1368 // create a random jet within the acceptance
1369 Double_t etaMax = fTrackEtaWindow - fRparam;
1372 Double_t pTC = 1; // small number
1373 for(int ir = 0;ir < fNRandomCones;ir++){
1374 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1375 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1377 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1378 Double_t pZC = pTC/TMath::Tan(thetaC);
1379 Double_t pXC = pTC * TMath::Cos(phiC);
1380 Double_t pYC = pTC * TMath::Sin(phiC);
1381 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1382 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
1384 for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1385 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1386 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1391 // test for overlap with previous cones to avoid double counting
1392 for(int iic = 0;iic<ir;iic++){
1393 AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1395 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1402 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1403 tmpRecC.SetPtLeading(-1.);
1404 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1405 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1406 }// loop over random cones creation
1409 // loop over the reconstructed particles and add up the pT in the random cones
1410 // maybe better to loop over randomized particles not in the real jets...
1411 // but this by definition brings dow average energy in the whole event
1412 AliAODJet vTmpRanR(1,0,0,1);
1413 for(int i = 0; i < recParticles.GetEntries(); i++){
1414 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1415 if(fTCARandomConesOut){
1416 for(int ir = 0;ir < fNRandomCones;ir++){
1417 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1418 if(jC&&jC->DeltaR(vp)<fRparam){
1419 if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1420 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1421 if(vp->Pt() > jC->GetPtLeading()) jC->SetPtLeading(vp->Pt());
1424 }// add up energy in cone
1426 // the randomized input changes eta and phi, but keeps the p_T
1427 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1428 Double_t pTR = vp->Pt();
1429 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1430 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1432 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1433 Double_t pZR = pTR/TMath::Tan(thetaR);
1435 Double_t pXR = pTR * TMath::Cos(phiR);
1436 Double_t pYR = pTR * TMath::Sin(phiR);
1437 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1438 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1439 if(fTCARandomConesOutRan){
1440 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1441 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1442 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1443 if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1444 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1445 if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt());
1450 }// loop over recparticles
1452 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1453 if(fTCARandomConesOut){
1454 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1455 // rescale the momentum vectors for the random cones
1457 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1459 Double_t etaC = rC->Eta();
1460 Double_t phiC = rC->Phi();
1461 // massless jet, unit vector
1462 pTC = rC->ChargedBgEnergy();
1463 if(pTC<=0)pTC = 0.001; // for almost empty events
1464 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1465 Double_t pZC = pTC/TMath::Tan(thetaC);
1466 Double_t pXC = pTC * TMath::Cos(phiC);
1467 Double_t pYC = pTC * TMath::Sin(phiC);
1468 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1469 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1470 rC->SetBgEnergy(0,0);
1471 rC->SetEffArea(jetArea,0);
1475 if(fTCARandomConesOutRan){
1476 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1477 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1480 Double_t etaC = rC->Eta();
1481 Double_t phiC = rC->Phi();
1482 // massless jet, unit vector
1483 pTC = rC->ChargedBgEnergy();
1484 if(pTC<=0)pTC = 0.001;// for almost empty events
1485 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1486 Double_t pZC = pTC/TMath::Tan(thetaC);
1487 Double_t pXC = pTC * TMath::Cos(phiC);
1488 Double_t pYC = pTC * TMath::Sin(phiC);
1489 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1490 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1491 rC->SetBgEnergy(0,0);
1492 rC->SetEffArea(jetArea,0);
1496 }// if(fNRandomCones
1498 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1500 if(fAODJetBackgroundOut){
1501 vector<fastjet::PseudoJet> jets2=sortedJets;
1502 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1506 Double_t meanarea1=0.;
1509 Double_t meanarea2=0.;
1513 Double_t meanarea4=0.;
1515 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1516 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1518 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1519 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1520 clustSeq.get_median_rho_and_sigma(sortedJets, range, true, bkg4, sigma4, meanarea4, true);
1521 fAODJetBackgroundOut->SetBackground(3,bkg4,sigma4,meanarea4);
1523 //////////////////////
1525 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1526 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1527 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1528 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1533 if(fStoreRhoLeadingTrackCorr) {
1534 vector<fastjet::PseudoJet> jets3=sortedJets;
1535 if(jets3.size()>2) jets3.erase(jets3.begin(),jets3.begin()+2);
1539 Double_t meanarea2=0.;
1541 clustSeq.get_median_rho_and_sigma(jets3, range, false, bkg2, sigma2, meanarea2, true);
1542 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1544 //Get leading track in event
1545 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1547 fh3CentvsRhoLeadingTrackPt->Fill(cent,bkg2,leading->Pt());
1548 fh3CentvsSigmaLeadingTrackPt->Fill(cent,sigma2,leading->Pt());
1549 fh3MultvsRhoLeadingTrackPt->Fill(nCh,bkg2,leading->Pt());
1550 fh3MultvsSigmaLeadingTrackPt->Fill(nCh,sigma2,leading->Pt());
1553 //Calculate rho with 4-vector area method for different quadrants with respect to the leading track in the event
1554 //exclude 2 leading jets from event
1555 //Quadrant 1: |phi_leadingTrack - phi_bkgCluster| < pi/2./2. - R (Near side to leading track)
1556 //Quadrant 2: pi/2 - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < pi/2 + (pi/2./2. - R) (Orthogonal to leading track)
1557 //Quadrant 3: pi - (pi/2/2 - R) < |phi_leadingTrack - phi_bkgCluster| < pi + (pi/2./2. - R) (Away side to leading track)
1558 //Quadrant 4: 3/2*pi - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < 3/2*pi + (pi/2./2. - R) (Orthogonal to leading track)
1560 //Assuming R=0.2 for background clusters
1562 Double_t bkg2Q[4] = {0.};
1563 Double_t sigma2Q[4] = {0.};
1564 Double_t meanarea2Q[4] = {0.};
1566 //Get phi of leading track in event
1567 Float_t phiLeadingTrack = leading->Phi();
1568 Float_t phiStep = (TMath::Pi()/2./2. - 0.2);
1570 //Quadrant1 - near side
1571 phiMin = phiLeadingTrack - phiStep;
1572 phiMax = phiLeadingTrack + phiStep;
1573 fastjet::RangeDefinition rangeQ1(rapMin,rapMax, phiMin, phiMax);
1574 clustSeq.get_median_rho_and_sigma(jets3, rangeQ1, false, bkg2Q[0], sigma2Q[0], meanarea2Q[0], true);
1576 //Quadrant2 - orthogonal
1577 Double_t phiQ2 = phiLeadingTrack + TMath::Pi()/2.;
1578 if(phiQ2>TMath::TwoPi()) phiQ2 = phiQ2 - TMath::TwoPi();
1579 phiMin = phiQ2 - phiStep;
1580 phiMax = phiQ2 + phiStep;
1581 fastjet::RangeDefinition rangeQ2(rapMin,rapMax, phiMin, phiMax);
1582 clustSeq.get_median_rho_and_sigma(jets3, rangeQ2, false, bkg2Q[1], sigma2Q[1], meanarea2Q[1], true);
1584 //Quadrant3 - away side
1585 Double_t phiQ3 = phiLeadingTrack + TMath::Pi();
1586 if(phiQ3>TMath::TwoPi()) phiQ3 = phiQ3 - TMath::TwoPi();
1587 phiMin = phiQ3 - phiStep;
1588 phiMax = phiQ3 + phiStep;
1589 fastjet::RangeDefinition rangeQ3(rapMin,rapMax, phiMin, phiMax);
1590 clustSeq.get_median_rho_and_sigma(jets3, rangeQ3, false, bkg2Q[2], sigma2Q[2], meanarea2Q[2], true);
1592 //Quadrant4 - othogonal
1593 Double_t phiQ4 = phiLeadingTrack + 3./2.*TMath::Pi();
1594 if(phiQ4>TMath::TwoPi()) phiQ4 = phiQ4 - TMath::TwoPi();
1595 phiMin = phiQ4 - phiStep;
1596 phiMax = phiQ4 + phiStep;
1597 fastjet::RangeDefinition rangeQ4(rapMin,rapMax, phiMin, phiMax);
1598 clustSeq.get_median_rho_and_sigma(jets3, rangeQ4, false, bkg2Q[3], sigma2Q[3], meanarea2Q[3], true);
1600 fh3CentvsRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0],leading->Pt());
1601 fh3CentvsRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1],leading->Pt());
1602 fh3CentvsRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2],leading->Pt());
1603 fh3CentvsRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3],leading->Pt());
1605 fh3CentvsSigmaLeadingTrackPtQ1->Fill(cent,sigma2Q[0],leading->Pt());
1606 fh3CentvsSigmaLeadingTrackPtQ2->Fill(cent,sigma2Q[1],leading->Pt());
1607 fh3CentvsSigmaLeadingTrackPtQ3->Fill(cent,sigma2Q[2],leading->Pt());
1608 fh3CentvsSigmaLeadingTrackPtQ4->Fill(cent,sigma2Q[3],leading->Pt());
1610 fh3MultvsRhoLeadingTrackPtQ1->Fill(nCh,bkg2Q[0],leading->Pt());
1611 fh3MultvsRhoLeadingTrackPtQ2->Fill(nCh,bkg2Q[1],leading->Pt());
1612 fh3MultvsRhoLeadingTrackPtQ3->Fill(nCh,bkg2Q[2],leading->Pt());
1613 fh3MultvsRhoLeadingTrackPtQ4->Fill(nCh,bkg2Q[3],leading->Pt());
1615 fh3MultvsSigmaLeadingTrackPtQ1->Fill(nCh,sigma2Q[0],leading->Pt());
1616 fh3MultvsSigmaLeadingTrackPtQ2->Fill(nCh,sigma2Q[1],leading->Pt());
1617 fh3MultvsSigmaLeadingTrackPtQ3->Fill(nCh,sigma2Q[2],leading->Pt());
1618 fh3MultvsSigmaLeadingTrackPtQ4->Fill(nCh,sigma2Q[3],leading->Pt());
1620 fh3CentvsDeltaRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0]-bkg2,leading->Pt());
1621 fh3CentvsDeltaRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1]-bkg2,leading->Pt());
1622 fh3CentvsDeltaRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2]-bkg2,leading->Pt());
1623 fh3CentvsDeltaRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3]-bkg2,leading->Pt());
1631 // fill track information
1632 Int_t nTrackOver = recParticles.GetSize();
1633 // do the same for tracks and jets
1636 TIterator *recIter = recParticles.MakeIterator();
1637 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1638 Float_t pt = tmpRec->Pt();
1640 // Printf("Leading track p_t %3.3E",pt);
1641 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1642 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1643 while(pt<ptCut&&tmpRec){
1645 tmpRec = (AliVParticle*)(recIter->Next());
1650 if(nTrackOver<=0)break;
1651 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1655 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1656 Float_t phi = leading->Phi();
1657 if(phi<0)phi+=TMath::Pi()*2.;
1658 Float_t eta = leading->Eta();
1660 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1661 Float_t tmpPt = tmpRec->Pt();
1662 Float_t tmpEta = tmpRec->Eta();
1663 fh1PtTracksRecIn->Fill(tmpPt);
1664 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1665 if(tmpRec==leading){
1666 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1667 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1671 Float_t tmpPhi = tmpRec->Phi();
1673 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1674 Float_t dPhi = phi - tmpPhi;
1675 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1676 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1677 Float_t dEta = eta - tmpRec->Eta();
1678 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1679 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1680 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1685 // find the random jets
1687 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1689 // fill the jet information from random track
1690 const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1691 const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1693 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1695 // loop over all jets an fill information, first on is the leading jet
1697 Int_t nRecOverRan = inclusiveJetsRan.size();
1698 Int_t nRecRan = inclusiveJetsRan.size();
1700 if(inclusiveJetsRan.size()>0){
1701 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1702 Float_t pt = leadingJet.Pt();
1705 TLorentzVector vecarearanb;
1707 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1708 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1709 while(pt<ptCut&&iCount<nRecRan){
1713 pt = sortedJetsRan[iCount].perp();
1716 if(nRecOverRan<=0)break;
1717 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1719 Float_t phi = leadingJet.Phi();
1720 if(phi<0)phi+=TMath::Pi()*2.;
1721 pt = leadingJet.Pt();
1723 // correlation of leading jet with random tracks
1725 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1727 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1729 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1730 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1731 Float_t dPhi = phi - tmpPhi;
1732 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1733 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1734 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1735 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1738 Int_t nAodOutJetsRan = 0;
1739 AliAODJet *aodOutJetRan = 0;
1740 for(int j = 0; j < nRecRan;j++){
1741 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1742 Float_t tmpPt = tmpRec.Pt();
1743 fh1PtJetsRecInRan->Fill(tmpPt);
1744 // Fill Spectra with constituents
1745 const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1746 fh1NConstRecRan->Fill(constituents.size());
1747 fh2NConstPtRan->Fill(tmpPt,constituents.size());
1748 fh2PtNchRan->Fill(nCh,tmpPt);
1749 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1752 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1753 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1754 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1755 aodOutJetRan->GetRefTracks()->Clear();
1756 aodOutJetRan->SetEffArea(arearan,0);
1757 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1758 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1759 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1764 Float_t tmpPhi = tmpRec.Phi();
1765 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1768 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1769 fh1NConstLeadingRecRan->Fill(constituents.size());
1770 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1776 if(fAODJetBackgroundOut){
1779 Double_t meanarea3=0.;
1780 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1781 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1782 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1784 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1785 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1794 // do the event selection if activated
1795 if(fJetTriggerPtCut>0){
1796 bool select = false;
1797 Float_t minPt = fJetTriggerPtCut;
1799 // hard coded for now ...
1800 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1801 if(cent<10)minPt = 50;
1802 else if(cent<30)minPt = 42;
1803 else if(cent<50)minPt = 28;
1804 else if(cent<80)minPt = 18;
1807 if(externalBackground)rho = externalBackground->GetBackground(2);
1809 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1810 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1811 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1820 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1821 fh1CentralitySelect->Fill(cent);
1822 fh1ZSelect->Fill(zVtx);
1823 aodH->SetFillAOD(kTRUE);
1827 if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1828 if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1829 if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1830 if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1832 PostData(1, fHistList);
1835 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1838 // Terminate analysis
1840 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1842 if(fMomResH1Fit) delete fMomResH1Fit;
1843 if(fMomResH2Fit) delete fMomResH2Fit;
1844 if(fMomResH3Fit) delete fMomResH3Fit;
1849 Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1852 // get list of tracks/particles for different types
1855 if(fDebug>2) Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1858 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly || type==kTrackAODMCextra || type==kTrackAODMCextraonly){
1860 if(type!=kTrackAODextraonly && type!=kTrackAODMCextraonly) {
1861 AliAODEvent *aod = 0;
1862 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1863 else aod = AODEvent();
1865 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1869 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1870 AliAODTrack *tr = aod->GetTrack(it);
1871 Bool_t bGood = false;
1872 if(fFilterType == 0)bGood = true;
1873 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1874 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1875 if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1876 if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1879 if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1880 if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1883 if(tr->Pt()<fTrackPtCut){
1884 if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1887 if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1892 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1893 AliAODEvent *aod = 0;
1894 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1895 else aod = AODEvent();
1900 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1901 if(!aodExtraTracks)return iCount;
1902 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1903 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1904 if (!track) continue;
1906 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1907 if(!trackAOD)continue;
1908 Bool_t bGood = false;
1909 if(fFilterType == 0)bGood = true;
1910 else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1911 else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1912 if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1913 if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1914 if(trackAOD->Pt()<fTrackPtCut) continue;
1915 if(fDebug) printf("pt extra track %.2f \n", trackAOD->Pt());
1916 list->Add(trackAOD);
1921 if(type==kTrackAODMCextra || type==kTrackAODMCextraonly) { //embed generator level particles
1922 AliAODEvent *aod = 0;
1923 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1924 else aod = AODEvent();
1928 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraMCparticles"));
1929 if(!aodExtraTracks)return iCount;
1930 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1931 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1933 if(fDebug) printf("track %d does not exist\n",it);
1938 mom[0] = track->Pt();
1939 mom[1] = track->Phi();
1940 mom[2] = track->Theta();
1942 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1944 if(fDebug) printf("trackAOD %d does not exist\n",it);
1948 trackAOD->SetFlags(AliESDtrack::kEmbedded);
1949 trackAOD->SetP(mom,kFALSE);
1950 trackAOD->SetFilterMap(fFilterMask);
1952 if(fDebug) printf("pt extra track %.2f \n", trackAOD->Pt());
1954 if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1955 if(trackAOD->Pt()<fTrackPtCut) continue;
1956 list->Add(trackAOD);
1962 else if (type == kTrackKineAll||type == kTrackKineCharged){
1963 AliMCEvent* mcEvent = MCEvent();
1964 if(!mcEvent)return iCount;
1965 // we want to have alivpartilces so use get track
1966 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1967 if(!mcEvent->IsPhysicalPrimary(it))continue;
1968 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1969 if(type == kTrackKineAll){
1970 if(part->Pt()<fTrackPtCut)continue;
1974 else if(type == kTrackKineCharged){
1975 if(part->Particle()->GetPDG()->Charge()==0)continue;
1976 if(part->Pt()<fTrackPtCut)continue;
1982 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1983 AliAODEvent *aod = 0;
1984 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1985 else aod = AODEvent();
1986 if(!aod)return iCount;
1987 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1988 if(!tca)return iCount;
1989 for(int it = 0;it < tca->GetEntriesFast();++it){
1990 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1991 if(!part->IsPhysicalPrimary())continue;
1992 if(type == kTrackAODMCAll){
1993 if(part->Pt()<fTrackPtCut)continue;
1997 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1998 if(part->Charge()==0)continue;
1999 if(part->Pt()<fTrackPtCut)continue;
2000 if(kTrackAODMCCharged){
2004 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
2015 void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() {
2018 AliInfo("Trying to connect to AliEn ...");
2019 TGrid::Connect("alien://");
2022 TFile *f = TFile::Open(fPathTrPtResolution.Data());
2026 TProfile *fProfPtPtSigma1PtGlobSt = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
2027 TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
2028 TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
2030 SetSmearResolution(kTRUE);
2031 SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD);
2036 void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() {
2039 AliInfo("Trying to connect to AliEn ...");
2040 TGrid::Connect("alien://");
2043 TFile *f = TFile::Open(fPathTrEfficiency.Data());
2046 TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
2047 TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
2048 TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
2050 SetDiceEfficiency(kTRUE);
2052 if(fChangeEfficiencyFraction>0.) {
2054 TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin");
2056 for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) {
2057 Double_t content = hEffPosGlobSt->GetBinContent(bin);
2058 hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction);
2061 SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
2065 SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
2069 void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
2072 // set mom res profiles
2075 if(fMomResH1) delete fMomResH1;
2076 if(fMomResH2) delete fMomResH2;
2077 if(fMomResH3) delete fMomResH3;
2079 fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
2080 fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
2081 fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
2085 void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
2087 // set tracking efficiency histos
2090 fhEffH1 = (TH1*)h1->Clone("fhEffH1");
2091 fhEffH2 = (TH1*)h2->Clone("fhEffH2");
2092 fhEffH3 = (TH1*)h3->Clone("fhEffH3");
2095 Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
2098 // Get smearing on generated momentum
2101 //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
2103 TProfile *fMomRes = 0x0;
2104 if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
2105 if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
2106 if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
2113 Double_t smear = 0.;
2116 if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
2117 if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
2118 if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
2122 Int_t bin = fMomRes->FindBin(pt);
2124 smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
2128 if(fMomRes) delete fMomRes;
2133 void AliAnalysisTaskJetCluster::FitMomentumResolution() {
2135 // Fit linear function on momentum resolution at high pT
2138 if(!fMomResH1Fit && fMomResH1) {
2139 fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
2140 fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
2141 fMomResH1Fit ->SetRange(5.,100.);
2144 if(!fMomResH2Fit && fMomResH2) {
2145 fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
2146 fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
2147 fMomResH2Fit ->SetRange(5.,100.);
2150 if(!fMomResH3Fit && fMomResH3) {
2151 fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
2152 fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
2153 fMomResH3Fit ->SetRange(5.,100.);
2159 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
2160 for(int i = 0; i < particles.GetEntries(); i++){
2161 AliVParticle *vp = (AliVParticle*)particles.At(i);
2162 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
2163 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
2164 jInp.set_user_index(i);
2165 inputParticles.push_back(jInp);