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),
111 fFilterMaskBestPt(0),
114 fTrackTypeRec(kTrackUndef),
115 fTrackTypeGen(kTrackUndef),
117 fNSkipLeadingCone(0),
121 fTrackEtaWindow(0.9),
123 fApplySharedClusterCut(0),
126 fJetOutputMinPt(0.150),
127 fMaxTrackPtInJet(100.),
133 fStoreRhoLeadingTrackCorr(kFALSE),
135 fBackgroundBranch(""),
146 fUseTrPtResolutionSmearing(kFALSE),
147 fUseDiceEfficiency(kFALSE),
148 fDiceEfficiencyMinPt(-1.),
149 fUseTrPtResolutionFromOADB(kFALSE),
150 fUseTrEfficiencyFromOADB(kFALSE),
151 fPathTrPtResolution(""),
152 fPathTrEfficiency(""),
153 fChangeEfficiencyFraction(0.),
154 fEfficiencyFixed(1.),
156 fAlgorithm(fastjet::kt_algorithm),
157 fStrategy(fastjet::Best),
158 fRecombScheme(fastjet::BIpt_scheme),
159 fAreaType(fastjet::active_area),
161 fActiveAreaRepeats(1),
165 fTCARandomConesOut(0x0),
166 fTCARandomConesOutRan(0x0),
167 fAODJetBackgroundOut(0x0),
173 fh1PtHardTrials(0x0),
176 fh1NConstLeadingRec(0x0),
178 fh1PtJetsLeadingRecIn(0x0),
179 fh1PtJetConstRec(0x0),
180 fh1PtJetConstLeadingRec(0x0),
181 fh1PtTracksRecIn(0x0),
182 fh1PtTracksLeadingRecIn(0x0),
184 fh1NConstRecRan(0x0),
185 fh1PtJetsLeadingRecInRan(0x0),
186 fh1NConstLeadingRecRan(0x0),
187 fh1PtJetsRecInRan(0x0),
188 fh1PtTracksGenIn(0x0),
190 fh1CentralityPhySel(0x0),
192 fh1CentralitySelect(0x0),
197 fh2NRecTracksPt(0x0),
199 fh2NConstLeadingPt(0x0),
201 fh2LeadingJetPhiEta(0x0),
203 fh2LeadingJetEtaPt(0x0),
205 fh2LeadingTrackEtaPt(0x0),
206 fh2JetsLeadingPhiEta(0x0),
207 fh2JetsLeadingPhiPt(0x0),
208 fh2TracksLeadingPhiEta(0x0),
209 fh2TracksLeadingPhiPt(0x0),
210 fh2TracksLeadingJetPhiPt(0x0),
211 fh2JetsLeadingPhiPtW(0x0),
212 fh2TracksLeadingPhiPtW(0x0),
213 fh2TracksLeadingJetPhiPtW(0x0),
214 fh2NRecJetsPtRan(0x0),
216 fh2NConstLeadingPtRan(0x0),
221 fh2TracksLeadingJetPhiPtRan(0x0),
222 fh2TracksLeadingJetPhiPtWRan(0x0),
223 fh3CentvsRhoLeadingTrackPt(0x0),
224 fh3CentvsSigmaLeadingTrackPt(0x0),
225 fh3MultvsRhoLeadingTrackPt(0x0),
226 fh3MultvsSigmaLeadingTrackPt(0x0),
227 fh3CentvsRhoLeadingTrackPtQ1(0x0),
228 fh3CentvsRhoLeadingTrackPtQ2(0x0),
229 fh3CentvsRhoLeadingTrackPtQ3(0x0),
230 fh3CentvsRhoLeadingTrackPtQ4(0x0),
231 fh3CentvsSigmaLeadingTrackPtQ1(0x0),
232 fh3CentvsSigmaLeadingTrackPtQ2(0x0),
233 fh3CentvsSigmaLeadingTrackPtQ3(0x0),
234 fh3CentvsSigmaLeadingTrackPtQ4(0x0),
235 fh3MultvsRhoLeadingTrackPtQ1(0x0),
236 fh3MultvsRhoLeadingTrackPtQ2(0x0),
237 fh3MultvsRhoLeadingTrackPtQ3(0x0),
238 fh3MultvsRhoLeadingTrackPtQ4(0x0),
239 fh3MultvsSigmaLeadingTrackPtQ1(0x0),
240 fh3MultvsSigmaLeadingTrackPtQ2(0x0),
241 fh3MultvsSigmaLeadingTrackPtQ3(0x0),
242 fh3MultvsSigmaLeadingTrackPtQ4(0x0),
243 fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0),
244 fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0),
245 fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0),
246 fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0),
247 fh2PtGenPtSmeared(0x0),
249 fp1PtResolution(0x0),
256 for(int i = 0;i<3;i++){
257 fh1BiARandomCones[i] = 0;
258 fh1BiARandomConesRan[i] = 0;
260 for(int i = 0;i<kMaxCent;i++){
261 fh2JetsLeadingPhiPtC[i] = 0;
262 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
263 fh2TracksLeadingJetPhiPtC[i] = 0;
264 fh2TracksLeadingJetPhiPtWC[i] = 0;
268 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
269 AliAnalysisTaskSE(name),
273 fUseAODTrackInput(kFALSE),
274 fUseAODMCInput(kFALSE),
275 fUseBackgroundCalc(kFALSE),
276 fEventSelection(kFALSE),
277 fRequireVZEROAC(kFALSE),
278 fRequireTZEROvtx(kFALSE),
281 fFilterMaskBestPt(0),
284 fTrackTypeRec(kTrackUndef),
285 fTrackTypeGen(kTrackUndef),
287 fNSkipLeadingCone(0),
291 fTrackEtaWindow(0.9),
293 fApplySharedClusterCut(0),
296 fJetOutputMinPt(0.150),
297 fMaxTrackPtInJet(100.),
303 fStoreRhoLeadingTrackCorr(kFALSE),
305 fBackgroundBranch(""),
316 fUseTrPtResolutionSmearing(kFALSE),
317 fUseDiceEfficiency(kFALSE),
318 fDiceEfficiencyMinPt(-1.),
319 fUseTrPtResolutionFromOADB(kFALSE),
320 fUseTrEfficiencyFromOADB(kFALSE),
321 fPathTrPtResolution(""),
322 fPathTrEfficiency(""),
323 fChangeEfficiencyFraction(0.),
324 fEfficiencyFixed(1.),
326 fAlgorithm(fastjet::kt_algorithm),
327 fStrategy(fastjet::Best),
328 fRecombScheme(fastjet::BIpt_scheme),
329 fAreaType(fastjet::active_area),
331 fActiveAreaRepeats(1),
335 fTCARandomConesOut(0x0),
336 fTCARandomConesOutRan(0x0),
337 fAODJetBackgroundOut(0x0),
343 fh1PtHardTrials(0x0),
346 fh1NConstLeadingRec(0x0),
348 fh1PtJetsLeadingRecIn(0x0),
349 fh1PtJetConstRec(0x0),
350 fh1PtJetConstLeadingRec(0x0),
351 fh1PtTracksRecIn(0x0),
352 fh1PtTracksLeadingRecIn(0x0),
354 fh1NConstRecRan(0x0),
355 fh1PtJetsLeadingRecInRan(0x0),
356 fh1NConstLeadingRecRan(0x0),
357 fh1PtJetsRecInRan(0x0),
358 fh1PtTracksGenIn(0x0),
360 fh1CentralityPhySel(0x0),
362 fh1CentralitySelect(0x0),
367 fh2NRecTracksPt(0x0),
369 fh2NConstLeadingPt(0x0),
371 fh2LeadingJetPhiEta(0x0),
373 fh2LeadingJetEtaPt(0x0),
375 fh2LeadingTrackEtaPt(0x0),
376 fh2JetsLeadingPhiEta(0x0),
377 fh2JetsLeadingPhiPt(0x0),
378 fh2TracksLeadingPhiEta(0x0),
379 fh2TracksLeadingPhiPt(0x0),
380 fh2TracksLeadingJetPhiPt(0x0),
381 fh2JetsLeadingPhiPtW(0x0),
382 fh2TracksLeadingPhiPtW(0x0),
383 fh2TracksLeadingJetPhiPtW(0x0),
384 fh2NRecJetsPtRan(0x0),
386 fh2NConstLeadingPtRan(0x0),
391 fh2TracksLeadingJetPhiPtRan(0x0),
392 fh2TracksLeadingJetPhiPtWRan(0x0),
393 fh3CentvsRhoLeadingTrackPt(0x0),
394 fh3CentvsSigmaLeadingTrackPt(0x0),
395 fh3MultvsRhoLeadingTrackPt(0x0),
396 fh3MultvsSigmaLeadingTrackPt(0x0),
397 fh3CentvsRhoLeadingTrackPtQ1(0x0),
398 fh3CentvsRhoLeadingTrackPtQ2(0x0),
399 fh3CentvsRhoLeadingTrackPtQ3(0x0),
400 fh3CentvsRhoLeadingTrackPtQ4(0x0),
401 fh3CentvsSigmaLeadingTrackPtQ1(0x0),
402 fh3CentvsSigmaLeadingTrackPtQ2(0x0),
403 fh3CentvsSigmaLeadingTrackPtQ3(0x0),
404 fh3CentvsSigmaLeadingTrackPtQ4(0x0),
405 fh3MultvsRhoLeadingTrackPtQ1(0x0),
406 fh3MultvsRhoLeadingTrackPtQ2(0x0),
407 fh3MultvsRhoLeadingTrackPtQ3(0x0),
408 fh3MultvsRhoLeadingTrackPtQ4(0x0),
409 fh3MultvsSigmaLeadingTrackPtQ1(0x0),
410 fh3MultvsSigmaLeadingTrackPtQ2(0x0),
411 fh3MultvsSigmaLeadingTrackPtQ3(0x0),
412 fh3MultvsSigmaLeadingTrackPtQ4(0x0),
413 fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0),
414 fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0),
415 fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0),
416 fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0),
417 fh2PtGenPtSmeared(0x0),
419 fp1PtResolution(0x0),
426 for(int i = 0;i<3;i++){
427 fh1BiARandomCones[i] = 0;
428 fh1BiARandomConesRan[i] = 0;
430 for(int i = 0;i<kMaxCent;i++){
431 fh2JetsLeadingPhiPtC[i] = 0;
432 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
433 fh2TracksLeadingJetPhiPtC[i] = 0;
434 fh2TracksLeadingJetPhiPtWC[i] = 0;
436 DefineOutput(1, TList::Class());
441 Bool_t AliAnalysisTaskJetCluster::Notify()
444 // Implemented Notify() to read the cross sections
445 // and number of trials from pyxsec.root
450 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
454 // Create the output container
457 fRandom = new TRandom3(0);
463 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
467 if(fNonStdBranch.Length()!=0)
469 // only create the output branch if we have a name
470 // Create a new branch for jets...
471 // -> cleared in the UserExec....
472 // here we can also have the case that the brnaches are written to a separate file
475 fTCAJetsOut = new TClonesArray("AliAODJet", 0);
476 fTCAJetsOut->SetName(fNonStdBranch.Data());
477 AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
480 if(fJetTypes&kJetRan){
481 fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
482 fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
483 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
484 if( fEfficiencyFixed < 1.)
485 fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dEffFixed%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.)));
487 fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dFr%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
489 AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
492 if(fUseBackgroundCalc){
493 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
494 fAODJetBackgroundOut = new AliAODJetEventBackground();
495 fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
496 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
497 if( fEfficiencyFixed < 1.)
498 fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dEffFixed%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.)));
500 fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dFr%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
502 AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
505 // create the branch for the random cones with the same R
506 TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
507 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
508 if( fEfficiencyFixed < 1.)
509 cName = Form("%sDetector%d%dEffFixed%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.),fNSkipLeadingCone);
511 cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone);
515 if(!AODEvent()->FindListObject(cName.Data())){
516 fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
517 fTCARandomConesOut->SetName(cName.Data());
518 AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
521 // create the branch with the random for the random cones on the random event
522 if(fJetTypes&kRCRan){
523 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
524 if(!AODEvent()->FindListObject(cName.Data())){
525 fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
526 fTCARandomConesOutRan->SetName(cName.Data());
527 AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
532 if(fNonStdFile.Length()!=0){
534 // case that we have an AOD extension we need to fetch the jets from the extended output
535 // we identify the extension aod event by looking for the branchname
536 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
537 // case that we have an AOD extension we need can fetch the background maybe from the extended output
538 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
543 if(!fHistList)fHistList = new TList();
544 fHistList->SetOwner();
545 PostData(1, fHistList); // post data in any case once
547 Bool_t oldStatus = TH1::AddDirectoryStatus();
548 TH1::AddDirectory(kFALSE);
553 const Int_t nBinPt = 100;
554 Double_t binLimitsPt[nBinPt+1];
555 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
557 binLimitsPt[iPt] = 0.0;
560 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
564 const Int_t nBinPhi = 90;
565 Double_t binLimitsPhi[nBinPhi+1];
566 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
568 binLimitsPhi[iPhi] = -1.*TMath::Pi();
571 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
577 const Int_t nBinEta = 40;
578 Double_t binLimitsEta[nBinEta+1];
579 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
581 binLimitsEta[iEta] = -2.0;
584 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
588 const int nChMax = 5000;
590 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
591 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
593 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
594 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
597 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
598 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
600 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
601 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
602 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
603 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
606 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
607 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
608 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
610 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
611 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
612 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
613 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
614 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
615 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
616 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
617 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
618 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
619 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
621 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
622 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
623 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
625 fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
626 fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
627 fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
629 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
630 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
631 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
635 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
636 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
637 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
638 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
640 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
641 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);
642 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);
643 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);
647 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
648 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
649 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
650 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
652 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
653 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
654 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
655 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
657 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
658 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
659 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
660 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
664 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
665 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
666 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
667 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
668 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
669 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
670 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
671 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
672 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
673 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
674 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
675 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
677 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
678 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
679 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
680 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
682 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
683 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
684 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
685 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
687 if(fStoreRhoLeadingTrackCorr) {
688 fh3CentvsRhoLeadingTrackPt = new TH3F("fh3CentvsRhoLeadingTrackPt","centrality vs background density full event; centrality; #rho", 50,0.,100.,500,0.,250.,100,0.,100.);
689 fh3CentvsSigmaLeadingTrackPt = new TH3F("fh3CentvsSigmaLeadingTrackPt","centrality vs sigma full event; centrality; #sigma(#rho)", 50,0.,100.,50,0.,50.,100,0.,100.);
690 fh3MultvsRhoLeadingTrackPt = new TH3F("fh3MultvsRhoLeadingTrackPt","multiplicity vs background density full event; multiplicity; #rho", 100,0.,5000.,500,0.,250.,100,0.,100.);
691 fh3MultvsSigmaLeadingTrackPt = new TH3F("fh3MultvsSigmaLeadingTrackPt","multiplicity vs sigma full event; multiplicity; #sigma(#rho)", 100,0.,5000.,50,0.,50.,100,0.,100.);
694 fh3CentvsRhoLeadingTrackPtQ1 = new TH3F("fh3CentvsRhoLeadingTrackPtQ1","centrality vs background density Q1; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
695 fh3CentvsRhoLeadingTrackPtQ2 = new TH3F("fh3CentvsRhoLeadingTrackPtQ2","centrality vs background density Q2; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
696 fh3CentvsRhoLeadingTrackPtQ3 = new TH3F("fh3CentvsRhoLeadingTrackPtQ3","centrality vs background density Q3; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
697 fh3CentvsRhoLeadingTrackPtQ4 = new TH3F("fh3CentvsRhoLeadingTrackPtQ4","centrality vs background density Q4; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
699 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.);
700 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.);
701 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.);
702 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.);
704 fh3MultvsRhoLeadingTrackPtQ1 = new TH3F("fh3MultvsRhoLeadingTrackPtQ1","multiplicity vs background density Q1; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
705 fh3MultvsRhoLeadingTrackPtQ2 = new TH3F("fh3MultvsRhoLeadingTrackPtQ2","multiplicity vs background density Q2; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
706 fh3MultvsRhoLeadingTrackPtQ3 = new TH3F("fh3MultvsRhoLeadingTrackPtQ3","multiplicity vs background density Q3; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
707 fh3MultvsRhoLeadingTrackPtQ4 = new TH3F("fh3MultvsRhoLeadingTrackPtQ4","multiplicity vs background density Q4; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
709 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.);
710 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.);
711 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.);
712 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.);
715 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.);
716 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.);
717 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.);
718 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.);
720 fHistList->Add(fh3CentvsRhoLeadingTrackPt);
721 fHistList->Add(fh3CentvsSigmaLeadingTrackPt);
722 fHistList->Add(fh3MultvsRhoLeadingTrackPt);
723 fHistList->Add(fh3MultvsSigmaLeadingTrackPt);
725 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ1);
726 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ2);
727 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ3);
728 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ4);
730 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ1);
731 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ2);
732 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ3);
733 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ4);
735 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ1);
736 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ2);
737 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ3);
738 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ4);
740 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ1);
741 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ2);
742 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ3);
743 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ4);
745 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ1);
746 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ2);
747 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ3);
748 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ4);
752 //Detector level effects histos
753 fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
755 fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
756 fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
758 fHistList->Add(fh2PtGenPtSmeared);
759 fHistList->Add(fp1Efficiency);
760 fHistList->Add(fp1PtResolution);
762 if(fNRandomCones>0&&fUseBackgroundCalc){
763 for(int i = 0;i<3;i++){
764 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
765 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
769 for(int i = 0;i < kMaxCent;i++){
770 fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
771 fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
772 fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
773 fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
776 const Int_t saveLevel = 3; // large save level more histos
778 fHistList->Add(fh1Xsec);
779 fHistList->Add(fh1Trials);
781 fHistList->Add(fh1NJetsRec);
782 fHistList->Add(fh1NConstRec);
783 fHistList->Add(fh1NConstLeadingRec);
784 fHistList->Add(fh1PtJetsRecIn);
785 fHistList->Add(fh1PtJetsLeadingRecIn);
786 fHistList->Add(fh1PtTracksRecIn);
787 fHistList->Add(fh1PtTracksLeadingRecIn);
788 fHistList->Add(fh1PtJetConstRec);
789 fHistList->Add(fh1PtJetConstLeadingRec);
790 fHistList->Add(fh1NJetsRecRan);
791 fHistList->Add(fh1NConstRecRan);
792 fHistList->Add(fh1PtJetsLeadingRecInRan);
793 fHistList->Add(fh1NConstLeadingRecRan);
794 fHistList->Add(fh1PtJetsRecInRan);
795 fHistList->Add(fh1Nch);
796 fHistList->Add(fh1Centrality);
797 fHistList->Add(fh1CentralitySelect);
798 fHistList->Add(fh1CentralityPhySel);
799 fHistList->Add(fh1Z);
800 fHistList->Add(fh1ZSelect);
801 fHistList->Add(fh1ZPhySel);
802 if(fNRandomCones>0&&fUseBackgroundCalc){
803 for(int i = 0;i<3;i++){
804 fHistList->Add(fh1BiARandomCones[i]);
805 fHistList->Add(fh1BiARandomConesRan[i]);
808 for(int i = 0;i < kMaxCent;i++){
809 fHistList->Add(fh2JetsLeadingPhiPtC[i]);
810 fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
811 fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
812 fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
815 fHistList->Add(fh2NRecJetsPt);
816 fHistList->Add(fh2NRecTracksPt);
817 fHistList->Add(fh2NConstPt);
818 fHistList->Add(fh2NConstLeadingPt);
819 fHistList->Add(fh2PtNch);
820 fHistList->Add(fh2PtNchRan);
821 fHistList->Add(fh2PtNchN);
822 fHistList->Add(fh2PtNchNRan);
823 fHistList->Add(fh2JetPhiEta);
824 fHistList->Add(fh2LeadingJetPhiEta);
825 fHistList->Add(fh2JetEtaPt);
826 fHistList->Add(fh2LeadingJetEtaPt);
827 fHistList->Add(fh2TrackEtaPt);
828 fHistList->Add(fh2LeadingTrackEtaPt);
829 fHistList->Add(fh2JetsLeadingPhiEta );
830 fHistList->Add(fh2JetsLeadingPhiPt);
831 fHistList->Add(fh2TracksLeadingPhiEta);
832 fHistList->Add(fh2TracksLeadingPhiPt);
833 fHistList->Add(fh2TracksLeadingJetPhiPt);
834 fHistList->Add(fh2JetsLeadingPhiPtW);
835 fHistList->Add(fh2TracksLeadingPhiPtW);
836 fHistList->Add(fh2TracksLeadingJetPhiPtW);
837 fHistList->Add(fh2NRecJetsPtRan);
838 fHistList->Add(fh2NConstPtRan);
839 fHistList->Add(fh2NConstLeadingPtRan);
840 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
841 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
844 // =========== Switch on Sumw2 for all histos ===========
845 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
846 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
851 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
854 TH1::AddDirectory(oldStatus);
857 void AliAnalysisTaskJetCluster::LocalInit()
863 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
865 if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB();
866 if(fUseTrEfficiencyFromOADB) LoadTrEfficiencyRootFileFromOADB();
869 FitMomentumResolution();
873 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
876 // handle and reset the output jet branch
878 if(fTCAJetsOut)fTCAJetsOut->Delete();
879 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
880 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
881 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
882 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
884 AliAODJetEventBackground* externalBackground = 0;
885 if(!externalBackground&&fBackgroundBranch.Length()){
886 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
887 if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
888 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
891 // Execute analysis for current event
893 AliESDEvent *fESD = 0;
894 if(fUseAODTrackInput){
895 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
897 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
903 // assume that the AOD is in the general output...
906 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
910 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
914 //Check if information is provided detector level effects
915 if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE;
916 if( fEfficiencyFixed < 1. ) {
917 if (!fUseDiceEfficiency)
918 fUseDiceEfficiency = 1; // 1 is the default; 2 can be set by user
921 if(!fhEffH1 || !fhEffH2 || !fhEffH3 ) fUseDiceEfficiency = kFALSE;
924 Bool_t selectEvent = false;
925 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
931 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
932 TString vtxTitle(vtxAOD->GetTitle());
933 zVtx = vtxAOD->GetZ();
935 cent = fAOD->GetHeader()->GetCentrality();
936 if(cent<10)cenClass = 0;
937 else if(cent<30)cenClass = 1;
938 else if(cent<50)cenClass = 2;
939 else if(cent<80)cenClass = 3;
940 if(physicsSelection){
941 fh1CentralityPhySel->Fill(cent);
942 fh1ZPhySel->Fill(zVtx);
946 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
947 Float_t yvtx = vtxAOD->GetY();
948 Float_t xvtx = vtxAOD->GetX();
949 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
950 if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
951 if(physicsSelection){
957 if(cent<fCentCutLo||cent>fCentCutUp){
969 const AliAODVZERO *vzero = fAOD->GetVZEROData();
971 if((vzero->GetTriggerChargeA()>0)&&(vzero->GetTriggerChargeC()>0)){
976 const AliAODTZERO *tzero = fAOD->GetTZEROData();
978 if(TMath::Abs(tzero->GetT0VertexRaw())<100){
983 if(fRequireVZEROAC&&fRequireTZEROvtx)selectEvent = selectEvent&&V0&&T0;
984 else if(fRequireTZEROvtx)selectEvent = selectEvent&&T0;
985 else if(fRequireVZEROAC)selectEvent = selectEvent&&V0;
989 PostData(1, fHistList);
992 fh1Centrality->Fill(cent);
994 fh1Trials->Fill("#sum{ntrials}",1);
997 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
999 // ==== General variables needed
1003 // we simply fetch the tracks/mc particles as a list of AliVParticles
1008 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
1009 Float_t nCh = recParticles.GetEntries();
1010 Float_t nGen=genParticles.GetEntries();
1012 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
1013 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
1014 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
1016 // find the jets....
1018 vector<fastjet::PseudoJet> inputParticlesRec;
1019 vector<fastjet::PseudoJet> inputParticlesRecRan;
1021 // Generate the random cones
1023 AliAODJet vTmpRan(1,0,0,1);
1024 for(int i = 0; i < recParticles.GetEntries(); i++){
1025 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1027 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1028 // we take total momentum here
1030 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
1031 //Add particles to fastjet in case we are not running toy model
1032 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
1033 jInp.set_user_index(i);
1034 inputParticlesRec.push_back(jInp);
1036 else if(fUseDiceEfficiency) {
1038 // Dice to decide if particle is kept or not - toy model for efficiency
1040 Double_t sumEff = 0.;
1042 Double_t eff[3] = {0.};
1045 Double_t rnd = fRandom->Uniform(1.);
1046 if( fEfficiencyFixed<1. ) {
1047 sumEff = fEfficiencyFixed;
1048 if (fUseDiceEfficiency == 2) {
1049 sumEff = (nCh - fEfficiencyFixed*nGen) / nCh; // rescale eff; fEfficiencyFixed is wrt to nGen, but dicing is fraction of nCh
1054 Double_t pTtmp = pT;
1055 if(pT>10.) pTtmp = 10.;
1056 Double_t eff1 = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
1057 Double_t eff2 = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
1058 Double_t eff3 = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));
1060 //Sort efficiencies from large to small
1061 if(eff1>eff2 && eff1>eff3) {
1076 else if(eff2>eff1 && eff2>eff3) {
1091 else if(eff3>eff1 && eff3>eff2) {
1107 sumEff = eff[0]+eff[1]+eff[2];
1109 fp1Efficiency->Fill(vp->Pt(),sumEff);
1110 if(rnd>sumEff && pT > fDiceEfficiencyMinPt) continue;
1112 if(fUseTrPtResolutionSmearing) {
1113 //Smear momentum of generated particle
1114 Double_t smear = 1.;
1115 //Select hybrid track category
1117 smear = GetMomentumSmearing(cat[2],pT);
1118 else if(rnd<=(eff[2]+eff[1]))
1119 smear = GetMomentumSmearing(cat[1],pT);
1121 smear = GetMomentumSmearing(cat[0],pT);
1123 fp1PtResolution->Fill(vp->Pt(),smear);
1125 Double_t sigma = vp->Pt()*smear;
1126 Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
1128 Double_t phi = vp->Phi();
1129 Double_t theta = 2.*TMath::ATan(TMath::Exp(-1.*(vp->Eta())));
1130 Double_t pX = pTrec * TMath::Cos(phi);
1131 Double_t pY = pTrec * TMath::Sin(phi);
1132 Double_t pZ = pTrec/TMath::Tan(theta);
1133 Double_t p=TMath::Sqrt(pTrec*pTrec+pZ*pZ);
1135 fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);
1137 fastjet::PseudoJet jInp(pX,pY,pZ,p);
1138 jInp.set_user_index(i);
1139 inputParticlesRec.push_back(jInp);
1143 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
1144 jInp.set_user_index(i);
1145 inputParticlesRec.push_back(jInp);
1151 // the randomized input changes eta and phi, but keeps the p_T
1152 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1153 Double_t pT = vp->Pt();
1154 Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
1155 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
1157 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
1158 Double_t pZ = pT/TMath::Tan(theta);
1160 Double_t pX = pT * TMath::Cos(phi);
1161 Double_t pY = pT * TMath::Sin(phi);
1162 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
1163 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
1165 jInpRan.set_user_index(i);
1166 inputParticlesRecRan.push_back(jInpRan);
1167 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
1170 // fill the tref array, only needed when we write out jets
1173 fRef->Delete(); // make sure to delete before placement new...
1174 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
1175 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
1178 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ...
1182 if(inputParticlesRec.size()==0){
1183 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
1184 PostData(1, fHistList);
1189 // employ setters for these...
1192 // now create the object that holds info about ghosts
1194 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
1195 // reduce CPU time...
1197 fActiveAreaRepeats = 0;
1201 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
1202 fastjet::AreaType areaType = fastjet::active_area;
1203 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
1204 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
1205 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
1207 //range where to compute background
1208 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
1210 phiMax = 2*TMath::Pi();
1211 rapMax = fGhostEtamax - fRparam;
1212 rapMin = - fGhostEtamax + fRparam;
1213 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
1216 const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
1217 const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
1220 fh1NJetsRec->Fill(sortedJets.size());
1222 // loop over all jets an fill information, first on is the leading jet
1224 Int_t nRecOver = inclusiveJets.size();
1225 Int_t nRec = inclusiveJets.size();
1226 if(inclusiveJets.size()>0){
1227 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
1228 Double_t area = clustSeq.area(sortedJets[0]);
1229 leadingJet.SetEffArea(area,0);
1230 Float_t pt = leadingJet.Pt();
1231 Int_t nAodOutJets = 0;
1232 Int_t nAodOutTracks = 0;
1233 AliAODJet *aodOutJet = 0;
1236 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
1237 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
1238 while(pt<ptCut&&iCount<nRec){
1242 pt = sortedJets[iCount].perp();
1245 if(nRecOver<=0)break;
1246 fh2NRecJetsPt->Fill(ptCut,nRecOver);
1248 Float_t phi = leadingJet.Phi();
1249 if(phi<0)phi+=TMath::Pi()*2.;
1250 Float_t eta = leadingJet.Eta();
1252 if(externalBackground){
1253 // carefull has to be filled in a task before
1254 // todo, ReArrange to the botom
1255 pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
1257 pt = leadingJet.Pt() - pTback;
1258 // correlation of leading jet with tracks
1259 TIterator *recIter = recParticles.MakeIterator();
1261 AliVParticle *tmpRecTrack = 0;
1262 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
1263 Float_t tmpPt = tmpRecTrack->Pt();
1265 Float_t tmpPhi = tmpRecTrack->Phi();
1266 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1267 Float_t dPhi = phi - tmpPhi;
1268 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1269 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1270 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
1271 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
1273 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
1274 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1280 TLorentzVector vecareab;
1281 for(int j = 0; j < nRec;j++){
1282 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
1285 Float_t tmpPt = tmpRec.Pt();
1287 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
1288 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
1289 aodOutJet->GetRefTracks()->Clear();
1290 Double_t area1 = clustSeq.area(sortedJets[j]);
1291 aodOutJet->SetEffArea(area1,0);
1292 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
1293 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
1294 aodOutJet->SetVectorAreaCharged(&vecareab);
1298 Float_t tmpPtBack = 0;
1299 if(externalBackground){
1300 // carefull has to be filled in a task before
1301 // todo, ReArrange to the botom
1302 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
1304 tmpPt = tmpPt - tmpPtBack;
1305 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
1307 fh1PtJetsRecIn->Fill(tmpPt);
1308 // Fill Spectra with constituentsemacs
1309 const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
1311 fh1NConstRec->Fill(constituents.size());
1312 fh2PtNch->Fill(nCh,tmpPt);
1313 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
1314 fh2NConstPt->Fill(tmpPt,constituents.size());
1315 // loop over constiutents and fill spectrum
1317 AliVParticle *partLead = 0;
1318 Float_t ptLead = -1;
1320 for(unsigned int ic = 0; ic < constituents.size();ic++){
1321 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1323 fh1PtJetConstRec->Fill(part->Pt());
1325 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1326 if(part->Pt()>fMaxTrackPtInJet){
1327 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1330 if(part->Pt()>ptLead){
1331 ptLead = part->Pt();
1334 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1337 AliAODTrack *aodT = 0;
1340 //set pT of leading constituent of jet
1341 aodOutJet->SetPtLeading(partLead->Pt());
1342 aodT = dynamic_cast<AliAODTrack*>(partLead);
1344 if(aodT->TestFilterBit(fFilterMaskBestPt)){
1345 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
1352 Float_t tmpPhi = tmpRec.Phi();
1353 Float_t tmpEta = tmpRec.Eta();
1354 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1356 fh1PtJetsLeadingRecIn->Fill(tmpPt);
1357 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1358 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1359 fh1NConstLeadingRec->Fill(constituents.size());
1360 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1363 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1364 fh2JetEtaPt->Fill(tmpEta,tmpPt);
1365 Float_t dPhi = phi - tmpPhi;
1366 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1367 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1368 Float_t dEta = eta - tmpRec.Eta();
1369 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1370 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1372 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1373 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1375 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1376 }// loop over reconstructed jets
1381 // Add the random cones...
1382 if(fNRandomCones>0&&fTCARandomConesOut){
1383 // create a random jet within the acceptance
1384 Double_t etaMax = fTrackEtaWindow - fRparam;
1387 Double_t pTC = 1; // small number
1388 for(int ir = 0;ir < fNRandomCones;ir++){
1389 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1390 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1392 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1393 Double_t pZC = pTC/TMath::Tan(thetaC);
1394 Double_t pXC = pTC * TMath::Cos(phiC);
1395 Double_t pYC = pTC * TMath::Sin(phiC);
1396 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1397 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
1399 for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1400 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1401 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1406 // test for overlap with previous cones to avoid double counting
1407 for(int iic = 0;iic<ir;iic++){
1408 AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1410 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1417 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1418 tmpRecC.SetPtLeading(-1.);
1419 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1420 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1421 }// loop over random cones creation
1424 // loop over the reconstructed particles and add up the pT in the random cones
1425 // maybe better to loop over randomized particles not in the real jets...
1426 // but this by definition brings dow average energy in the whole event
1427 AliAODJet vTmpRanR(1,0,0,1);
1428 for(int i = 0; i < recParticles.GetEntries(); i++){
1429 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1430 if(fTCARandomConesOut){
1431 for(int ir = 0;ir < fNRandomCones;ir++){
1432 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1433 if(jC&&jC->DeltaR(vp)<fRparam){
1434 if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1435 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1436 if(vp->Pt() > jC->GetPtLeading()) jC->SetPtLeading(vp->Pt());
1439 }// add up energy in cone
1441 // the randomized input changes eta and phi, but keeps the p_T
1442 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1443 Double_t pTR = vp->Pt();
1444 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1445 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1447 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1448 Double_t pZR = pTR/TMath::Tan(thetaR);
1450 Double_t pXR = pTR * TMath::Cos(phiR);
1451 Double_t pYR = pTR * TMath::Sin(phiR);
1452 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1453 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1454 if(fTCARandomConesOutRan){
1455 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1456 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1457 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1458 if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1459 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1460 if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt());
1465 }// loop over recparticles
1467 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1468 if(fTCARandomConesOut){
1469 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1470 // rescale the momentum vectors for the random cones
1472 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1474 Double_t etaC = rC->Eta();
1475 Double_t phiC = rC->Phi();
1476 // massless jet, unit vector
1477 pTC = rC->ChargedBgEnergy();
1478 if(pTC<=0)pTC = 0.001; // for almost empty events
1479 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1480 Double_t pZC = pTC/TMath::Tan(thetaC);
1481 Double_t pXC = pTC * TMath::Cos(phiC);
1482 Double_t pYC = pTC * TMath::Sin(phiC);
1483 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1484 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1485 rC->SetBgEnergy(0,0);
1486 rC->SetEffArea(jetArea,0);
1490 if(fTCARandomConesOutRan){
1491 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1492 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1495 Double_t etaC = rC->Eta();
1496 Double_t phiC = rC->Phi();
1497 // massless jet, unit vector
1498 pTC = rC->ChargedBgEnergy();
1499 if(pTC<=0)pTC = 0.001;// for almost empty events
1500 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1501 Double_t pZC = pTC/TMath::Tan(thetaC);
1502 Double_t pXC = pTC * TMath::Cos(phiC);
1503 Double_t pYC = pTC * TMath::Sin(phiC);
1504 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1505 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1506 rC->SetBgEnergy(0,0);
1507 rC->SetEffArea(jetArea,0);
1511 }// if(fNRandomCones
1513 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1515 if(fAODJetBackgroundOut){
1516 vector<fastjet::PseudoJet> jets2=sortedJets;
1517 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1521 Double_t meanarea1=0.;
1524 Double_t meanarea2=0.;
1528 Double_t meanarea4=0.;
1530 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1531 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1533 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1534 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1535 clustSeq.get_median_rho_and_sigma(sortedJets, range, true, bkg4, sigma4, meanarea4, true);
1536 fAODJetBackgroundOut->SetBackground(3,bkg4,sigma4,meanarea4);
1538 //////////////////////
1540 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1541 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1542 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1543 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1548 if(fStoreRhoLeadingTrackCorr) {
1549 vector<fastjet::PseudoJet> jets3=sortedJets;
1550 if(jets3.size()>2) jets3.erase(jets3.begin(),jets3.begin()+2);
1554 Double_t meanarea2=0.;
1556 clustSeq.get_median_rho_and_sigma(jets3, range, false, bkg2, sigma2, meanarea2, true);
1557 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1559 //Get leading track in event
1560 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1562 fh3CentvsRhoLeadingTrackPt->Fill(cent,bkg2,leading->Pt());
1563 fh3CentvsSigmaLeadingTrackPt->Fill(cent,sigma2,leading->Pt());
1564 fh3MultvsRhoLeadingTrackPt->Fill(nCh,bkg2,leading->Pt());
1565 fh3MultvsSigmaLeadingTrackPt->Fill(nCh,sigma2,leading->Pt());
1568 //Calculate rho with 4-vector area method for different quadrants with respect to the leading track in the event
1569 //exclude 2 leading jets from event
1570 //Quadrant 1: |phi_leadingTrack - phi_bkgCluster| < pi/2./2. - R (Near side to leading track)
1571 //Quadrant 2: pi/2 - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < pi/2 + (pi/2./2. - R) (Orthogonal to leading track)
1572 //Quadrant 3: pi - (pi/2/2 - R) < |phi_leadingTrack - phi_bkgCluster| < pi + (pi/2./2. - R) (Away side to leading track)
1573 //Quadrant 4: 3/2*pi - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < 3/2*pi + (pi/2./2. - R) (Orthogonal to leading track)
1575 //Assuming R=0.2 for background clusters
1577 Double_t bkg2Q[4] = {0.};
1578 Double_t sigma2Q[4] = {0.};
1579 Double_t meanarea2Q[4] = {0.};
1581 //Get phi of leading track in event
1582 Float_t phiLeadingTrack = leading->Phi();
1583 Float_t phiStep = (TMath::Pi()/2./2. - 0.2);
1585 //Quadrant1 - near side
1586 phiMin = phiLeadingTrack - phiStep;
1587 phiMax = phiLeadingTrack + phiStep;
1588 fastjet::RangeDefinition rangeQ1(rapMin,rapMax, phiMin, phiMax);
1589 clustSeq.get_median_rho_and_sigma(jets3, rangeQ1, false, bkg2Q[0], sigma2Q[0], meanarea2Q[0], true);
1591 //Quadrant2 - orthogonal
1592 Double_t phiQ2 = phiLeadingTrack + TMath::Pi()/2.;
1593 if(phiQ2>TMath::TwoPi()) phiQ2 = phiQ2 - TMath::TwoPi();
1594 phiMin = phiQ2 - phiStep;
1595 phiMax = phiQ2 + phiStep;
1596 fastjet::RangeDefinition rangeQ2(rapMin,rapMax, phiMin, phiMax);
1597 clustSeq.get_median_rho_and_sigma(jets3, rangeQ2, false, bkg2Q[1], sigma2Q[1], meanarea2Q[1], true);
1599 //Quadrant3 - away side
1600 Double_t phiQ3 = phiLeadingTrack + TMath::Pi();
1601 if(phiQ3>TMath::TwoPi()) phiQ3 = phiQ3 - TMath::TwoPi();
1602 phiMin = phiQ3 - phiStep;
1603 phiMax = phiQ3 + phiStep;
1604 fastjet::RangeDefinition rangeQ3(rapMin,rapMax, phiMin, phiMax);
1605 clustSeq.get_median_rho_and_sigma(jets3, rangeQ3, false, bkg2Q[2], sigma2Q[2], meanarea2Q[2], true);
1607 //Quadrant4 - othogonal
1608 Double_t phiQ4 = phiLeadingTrack + 3./2.*TMath::Pi();
1609 if(phiQ4>TMath::TwoPi()) phiQ4 = phiQ4 - TMath::TwoPi();
1610 phiMin = phiQ4 - phiStep;
1611 phiMax = phiQ4 + phiStep;
1612 fastjet::RangeDefinition rangeQ4(rapMin,rapMax, phiMin, phiMax);
1613 clustSeq.get_median_rho_and_sigma(jets3, rangeQ4, false, bkg2Q[3], sigma2Q[3], meanarea2Q[3], true);
1615 fh3CentvsRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0],leading->Pt());
1616 fh3CentvsRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1],leading->Pt());
1617 fh3CentvsRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2],leading->Pt());
1618 fh3CentvsRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3],leading->Pt());
1620 fh3CentvsSigmaLeadingTrackPtQ1->Fill(cent,sigma2Q[0],leading->Pt());
1621 fh3CentvsSigmaLeadingTrackPtQ2->Fill(cent,sigma2Q[1],leading->Pt());
1622 fh3CentvsSigmaLeadingTrackPtQ3->Fill(cent,sigma2Q[2],leading->Pt());
1623 fh3CentvsSigmaLeadingTrackPtQ4->Fill(cent,sigma2Q[3],leading->Pt());
1625 fh3MultvsRhoLeadingTrackPtQ1->Fill(nCh,bkg2Q[0],leading->Pt());
1626 fh3MultvsRhoLeadingTrackPtQ2->Fill(nCh,bkg2Q[1],leading->Pt());
1627 fh3MultvsRhoLeadingTrackPtQ3->Fill(nCh,bkg2Q[2],leading->Pt());
1628 fh3MultvsRhoLeadingTrackPtQ4->Fill(nCh,bkg2Q[3],leading->Pt());
1630 fh3MultvsSigmaLeadingTrackPtQ1->Fill(nCh,sigma2Q[0],leading->Pt());
1631 fh3MultvsSigmaLeadingTrackPtQ2->Fill(nCh,sigma2Q[1],leading->Pt());
1632 fh3MultvsSigmaLeadingTrackPtQ3->Fill(nCh,sigma2Q[2],leading->Pt());
1633 fh3MultvsSigmaLeadingTrackPtQ4->Fill(nCh,sigma2Q[3],leading->Pt());
1635 fh3CentvsDeltaRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0]-bkg2,leading->Pt());
1636 fh3CentvsDeltaRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1]-bkg2,leading->Pt());
1637 fh3CentvsDeltaRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2]-bkg2,leading->Pt());
1638 fh3CentvsDeltaRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3]-bkg2,leading->Pt());
1646 // fill track information
1647 Int_t nTrackOver = recParticles.GetSize();
1648 // do the same for tracks and jets
1651 TIterator *recIter = recParticles.MakeIterator();
1652 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1653 Float_t pt = tmpRec->Pt();
1655 // Printf("Leading track p_t %3.3E",pt);
1656 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1657 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1658 while(pt<ptCut&&tmpRec){
1660 tmpRec = (AliVParticle*)(recIter->Next());
1665 if(nTrackOver<=0)break;
1666 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1670 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1671 Float_t phi = leading->Phi();
1672 if(phi<0)phi+=TMath::Pi()*2.;
1673 Float_t eta = leading->Eta();
1675 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1676 Float_t tmpPt = tmpRec->Pt();
1677 Float_t tmpEta = tmpRec->Eta();
1678 fh1PtTracksRecIn->Fill(tmpPt);
1679 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1680 if(tmpRec==leading){
1681 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1682 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1686 Float_t tmpPhi = tmpRec->Phi();
1688 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1689 Float_t dPhi = phi - tmpPhi;
1690 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1691 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1692 Float_t dEta = eta - tmpRec->Eta();
1693 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1694 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1695 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1700 // find the random jets
1702 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1704 // fill the jet information from random track
1705 const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1706 const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1708 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1710 // loop over all jets an fill information, first on is the leading jet
1712 Int_t nRecOverRan = inclusiveJetsRan.size();
1713 Int_t nRecRan = inclusiveJetsRan.size();
1715 if(inclusiveJetsRan.size()>0){
1716 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1717 Float_t pt = leadingJet.Pt();
1720 TLorentzVector vecarearanb;
1722 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1723 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1724 while(pt<ptCut&&iCount<nRecRan){
1728 pt = sortedJetsRan[iCount].perp();
1731 if(nRecOverRan<=0)break;
1732 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1734 Float_t phi = leadingJet.Phi();
1735 if(phi<0)phi+=TMath::Pi()*2.;
1736 pt = leadingJet.Pt();
1738 // correlation of leading jet with random tracks
1740 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1742 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1744 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1745 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1746 Float_t dPhi = phi - tmpPhi;
1747 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1748 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1749 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1750 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1753 Int_t nAodOutJetsRan = 0;
1754 AliAODJet *aodOutJetRan = 0;
1755 for(int j = 0; j < nRecRan;j++){
1756 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1757 Float_t tmpPt = tmpRec.Pt();
1758 fh1PtJetsRecInRan->Fill(tmpPt);
1759 // Fill Spectra with constituents
1760 const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1761 fh1NConstRecRan->Fill(constituents.size());
1762 fh2NConstPtRan->Fill(tmpPt,constituents.size());
1763 fh2PtNchRan->Fill(nCh,tmpPt);
1764 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1767 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1768 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1769 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1770 aodOutJetRan->GetRefTracks()->Clear();
1771 aodOutJetRan->SetEffArea(arearan,0);
1772 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1773 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1774 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1779 Float_t tmpPhi = tmpRec.Phi();
1780 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1783 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1784 fh1NConstLeadingRecRan->Fill(constituents.size());
1785 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1791 if(fAODJetBackgroundOut){
1794 Double_t meanarea3=0.;
1795 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1796 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1797 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1799 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1800 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1809 // do the event selection if activated
1810 if(fJetTriggerPtCut>0){
1811 bool select = false;
1812 Float_t minPt = fJetTriggerPtCut;
1814 // hard coded for now ...
1815 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1816 if(cent<10)minPt = 50;
1817 else if(cent<30)minPt = 42;
1818 else if(cent<50)minPt = 28;
1819 else if(cent<80)minPt = 18;
1822 if(externalBackground)rho = externalBackground->GetBackground(2);
1824 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1825 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1826 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1835 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1836 fh1CentralitySelect->Fill(cent);
1837 fh1ZSelect->Fill(zVtx);
1838 aodH->SetFillAOD(kTRUE);
1842 if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1843 if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1844 if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1845 if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1847 PostData(1, fHistList);
1850 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1853 // Terminate analysis
1855 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1857 if(fMomResH1Fit) delete fMomResH1Fit;
1858 if(fMomResH2Fit) delete fMomResH2Fit;
1859 if(fMomResH3Fit) delete fMomResH3Fit;
1864 Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1867 // get list of tracks/particles for different types
1870 if(fDebug>2) Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1873 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly || type==kTrackAODMCextra || type==kTrackAODMCextraonly){
1875 if(type!=kTrackAODextraonly && type!=kTrackAODMCextraonly) {
1876 AliAODEvent *aod = 0;
1877 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1878 else aod = AODEvent();
1880 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1884 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1885 AliAODTrack *tr = aod->GetTrack(it);
1886 Bool_t bGood = false;
1887 if(fFilterType == 0)bGood = true;
1888 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1889 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1890 if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1891 if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1895 // heavy flavor jets
1896 if(fFilterMask==528 && fUseHFcuts){
1897 Double_t ntpcClus = tr->GetTPCNcls();
1898 Double_t trPt=tr->Pt();
1899 TFormula NTPCClsCut("f1NClustersTPCLinearPtDep","70.+30./20.*x");
1901 if (trPt <= 20. && (ntpcClus < NTPCClsCut.Eval(trPt))) continue;
1902 else if (trPt > 20. && ntpcClus < 100) continue;
1904 if (AvoidDoubleCountingHF(aod,tr)) continue;
1908 if(fRequireITSRefit){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
1909 if (fApplySharedClusterCut) {
1910 Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
1911 if (frac > 0.4) continue;
1913 if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1914 if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1917 if(tr->Pt()<fTrackPtCut){
1918 if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1921 if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1926 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1927 AliAODEvent *aod = 0;
1928 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1929 else aod = AODEvent();
1934 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1935 if(!aodExtraTracks)return iCount;
1936 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1937 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1938 if (!track) continue;
1940 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1941 if(!trackAOD)continue;
1942 Bool_t bGood = false;
1943 if(fFilterType == 0)bGood = true;
1944 else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1945 else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1946 if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1947 if(fRequireITSRefit){if((trackAOD->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
1948 if (fApplySharedClusterCut) {
1949 Double_t frac = Double_t(trackAOD->GetTPCnclsS()) /Double_t(trackAOD->GetTPCncls());
1950 if (frac > 0.4) continue;
1954 if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1955 if(trackAOD->Pt()<fTrackPtCut) continue;
1956 if(fDebug) printf("pt extra track %.2f \n", trackAOD->Pt());
1957 list->Add(trackAOD);
1962 if(type==kTrackAODMCextra || type==kTrackAODMCextraonly) { //embed generator level particles
1963 AliAODEvent *aod = 0;
1964 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1965 else aod = AODEvent();
1969 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraMCparticles"));
1970 if(!aodExtraTracks)return iCount;
1971 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1972 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1973 AliAODMCParticle *partmc = dynamic_cast<AliAODMCParticle*> ((*aodExtraTracks)[it]);
1975 if(fDebug>10) printf("track %d does not exist\n",it);
1979 if(!partmc) continue;
1980 if(!partmc->IsPhysicalPrimary())continue;
1982 if (track->Pt()<fTrackPtCut) {
1983 if(fDebug>10) printf("track %d has too low pt %.2f\n",it,track->Pt());
1988 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*>((*aodExtraTracks)[it]);//(track);
1991 if(fDebug>10) printf("trackAOD %d does not exist\n",it);
1995 trackAOD->SetFlags(AliESDtrack::kEmbedded);
1996 trackAOD->SetFilterMap(fFilterMask);
1998 if(fDebug>10) printf("pt extra track %.2f \n", track->Pt());
2000 if(TMath::Abs(track->Eta())>fTrackEtaWindow) continue;
2001 if(track->Pt()<fTrackPtCut) continue;
2009 else if (type == kTrackKineAll||type == kTrackKineCharged){
2010 AliMCEvent* mcEvent = MCEvent();
2011 if(!mcEvent)return iCount;
2012 // we want to have alivpartilces so use get track
2013 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
2014 if(!mcEvent->IsPhysicalPrimary(it))continue;
2015 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
2016 if(type == kTrackKineAll){
2017 if(part->Pt()<fTrackPtCut)continue;
2021 else if(type == kTrackKineCharged){
2022 if(part->Particle()->GetPDG()->Charge()==0)continue;
2023 if(part->Pt()<fTrackPtCut)continue;
2029 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
2030 AliAODEvent *aod = 0;
2031 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
2032 else aod = AODEvent();
2033 if(!aod)return iCount;
2034 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
2035 if(!tca)return iCount;
2036 for(int it = 0;it < tca->GetEntriesFast();++it){
2037 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
2038 if(!part->IsPhysicalPrimary())continue;
2039 if(type == kTrackAODMCAll){
2040 if(part->Pt()<fTrackPtCut)continue;
2044 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
2045 if(part->Charge()==0)continue;
2046 if(part->Pt()<fTrackPtCut)continue;
2047 if(kTrackAODMCCharged){
2051 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
2058 else if (type == kTrackAODMCHF){
2060 AliAODEvent *aod = 0;
2061 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
2062 else aod = AODEvent();
2063 if(!aod)return iCount;
2064 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
2065 if(!tca)return iCount;
2066 for(int it = 0;it < tca->GetEntriesFast();++it){
2067 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
2069 int partpdg= part->PdgCode();
2070 if(!part->IsPhysicalPrimary() && !IsBMeson(partpdg) && !IsDMeson(partpdg) )continue;
2072 if (IsBMeson(partpdg) || IsDMeson(partpdg)) {
2073 iCount+= AddDaughters( list , part,tca);
2077 if(part->Pt()<fTrackPtCut) continue;
2078 if(TMath::Abs(part->Eta())>fTrackEtaWindow) continue;
2079 if(part->Charge()==0) continue;
2081 if((part->Pt()>=fTrackPtCut) && (TMath::Abs(part->Eta())<=fTrackEtaWindow) && (part->Charge()!=0))list->Add(part);
2091 Int_t AliAnalysisTaskJetCluster::AddDaughters(TList * list, AliAODMCParticle *part, TClonesArray * tca){
2093 Int_t nDaugthers = part->GetNDaughters();
2094 for(Int_t i=0;i< nDaugthers;++i){
2095 AliAODMCParticle *partdaughter = (AliAODMCParticle*)(tca->At(i));
2096 if(!partdaughter) continue;
2097 if(partdaughter->Pt()<fTrackPtCut)continue;
2098 if(TMath::Abs(partdaughter->Eta())>fTrackEtaWindow)continue;
2099 if(partdaughter->Charge()==0)continue;
2101 if(!IsDMeson(partdaughter->PdgCode()) && !IsBMeson(partdaughter->PdgCode()) ){
2102 list->Add(partdaughter);
2105 else count+=AddDaughters(list,part,tca);
2111 Bool_t AliAnalysisTaskJetCluster::AvoidDoubleCountingHF(AliAODEvent *aod, AliAODTrack *tr1){
2113 if(!(tr1->TestFilterBit(BIT(9)))) return kFALSE;
2115 Int_t idtr1 = tr1->GetID();
2117 for(int jt = 0;jt < aod->GetNumberOfTracks();++jt){
2119 const AliAODTrack *tr2 = aod->GetTrack(jt);
2120 Int_t idtr2 = tr2->GetID();
2122 if (!(tr2->TestFilterBit(BIT(4)))) continue;
2123 if (idtr1==(idtr2+1)*-1.) return kTRUE;
2128 void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() {
2131 AliInfo("Trying to connect to AliEn ...");
2132 TGrid::Connect("alien://");
2135 TFile *f = TFile::Open(fPathTrPtResolution.Data());
2139 TProfile *fProfPtPtSigma1PtGlobSt = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
2140 TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
2141 TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
2143 SetSmearResolution(kTRUE);
2144 SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD);
2149 void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() {
2152 AliInfo("Trying to connect to AliEn ...");
2153 TGrid::Connect("alien://");
2156 TFile *f = TFile::Open(fPathTrEfficiency.Data());
2159 TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
2160 TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
2161 TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
2163 SetDiceEfficiency(kTRUE);
2165 if(fChangeEfficiencyFraction>0.) {
2167 TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin");
2169 for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) {
2170 Double_t content = hEffPosGlobSt->GetBinContent(bin);
2171 hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction);
2174 SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
2178 SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
2182 void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
2185 // set mom res profiles
2188 if(fMomResH1) delete fMomResH1;
2189 if(fMomResH2) delete fMomResH2;
2190 if(fMomResH3) delete fMomResH3;
2192 fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
2193 fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
2194 fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
2198 void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
2200 // set tracking efficiency histos
2203 fhEffH1 = (TH1*)h1->Clone("fhEffH1");
2204 fhEffH2 = (TH1*)h2->Clone("fhEffH2");
2205 fhEffH3 = (TH1*)h3->Clone("fhEffH3");
2208 Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
2211 // Get smearing on generated momentum
2214 //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
2216 TProfile *fMomRes = 0x0;
2217 if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
2218 if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
2219 if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
2226 Double_t smear = 0.;
2229 if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
2230 if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
2231 if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
2235 Int_t bin = fMomRes->FindBin(pt);
2237 smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
2241 if(fMomRes) delete fMomRes;
2246 void AliAnalysisTaskJetCluster::FitMomentumResolution() {
2248 // Fit linear function on momentum resolution at high pT
2251 if(!fMomResH1Fit && fMomResH1) {
2252 fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
2253 fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
2254 fMomResH1Fit ->SetRange(5.,100.);
2257 if(!fMomResH2Fit && fMomResH2) {
2258 fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
2259 fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
2260 fMomResH2Fit ->SetRange(5.,100.);
2263 if(!fMomResH3Fit && fMomResH3) {
2264 fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
2265 fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
2266 fMomResH3Fit ->SetRange(5.,100.);
2272 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
2273 for(int i = 0; i < particles.GetEntries(); i++){
2274 AliVParticle *vp = (AliVParticle*)particles.At(i);
2275 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
2276 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
2277 jInp.set_user_index(i);
2278 inputParticles.push_back(jInp);
2287 bool AliAnalysisTaskJetCluster::IsBMeson(int pc){
2288 int bPdG[] = {511,521,10511,10521,513,523,10513,10523,20513,20523,20513,20523,515,525,531,
2289 10531,533,10533,20533,535,541,10541,543,10543,20543,545,551,10551,100551,
2290 110551,200551,210551,553,10553,20553,30553,100553,110553,120553,130553,200553,210553,220553,
2291 300553,9000533,9010553,555,10555,20555,100555,110555,120555,200555,557,100557};
2292 for(int i=0;i< (int)(sizeof(bPdG)/sizeof(int));++i) if(abs(pc) == bPdG[i]) return true;
2295 bool AliAnalysisTaskJetCluster::IsDMeson(int pc){
2296 int bPdG[] = {411,421,10411,10421,413,423,10413,10423,20431,20423,415,
2297 425,431,10431,433,10433,20433,435,441,10441,100441,443,10443,20443,
2298 100443,30443,9000443,9010443,9020443,445,100445};
2299 for(int i=0;i< (int)(sizeof(bPdG)/sizeof(int));++i) if(abs(pc) == bPdG[i]) return true;