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),
122 fApplySharedClusterCut(0),
125 fJetOutputMinPt(0.150),
126 fMaxTrackPtInJet(100.),
132 fStoreRhoLeadingTrackCorr(kFALSE),
134 fBackgroundBranch(""),
145 fUseTrPtResolutionSmearing(kFALSE),
146 fUseDiceEfficiency(kFALSE),
147 fDiceEfficiencyMinPt(-1.),
148 fUseTrPtResolutionFromOADB(kFALSE),
149 fUseTrEfficiencyFromOADB(kFALSE),
150 fPathTrPtResolution(""),
151 fPathTrEfficiency(""),
152 fChangeEfficiencyFraction(0.),
153 fEfficiencyFixed(1.),
155 fAlgorithm(fastjet::kt_algorithm),
156 fStrategy(fastjet::Best),
157 fRecombScheme(fastjet::BIpt_scheme),
158 fAreaType(fastjet::active_area),
160 fActiveAreaRepeats(1),
164 fTCARandomConesOut(0x0),
165 fTCARandomConesOutRan(0x0),
166 fAODJetBackgroundOut(0x0),
172 fh1PtHardTrials(0x0),
175 fh1NConstLeadingRec(0x0),
177 fh1PtJetsLeadingRecIn(0x0),
178 fh1PtJetConstRec(0x0),
179 fh1PtJetConstLeadingRec(0x0),
180 fh1PtTracksRecIn(0x0),
181 fh1PtTracksLeadingRecIn(0x0),
183 fh1NConstRecRan(0x0),
184 fh1PtJetsLeadingRecInRan(0x0),
185 fh1NConstLeadingRecRan(0x0),
186 fh1PtJetsRecInRan(0x0),
187 fh1PtTracksGenIn(0x0),
189 fh1CentralityPhySel(0x0),
191 fh1CentralitySelect(0x0),
196 fh2NRecTracksPt(0x0),
198 fh2NConstLeadingPt(0x0),
200 fh2LeadingJetPhiEta(0x0),
202 fh2LeadingJetEtaPt(0x0),
204 fh2LeadingTrackEtaPt(0x0),
205 fh2JetsLeadingPhiEta(0x0),
206 fh2JetsLeadingPhiPt(0x0),
207 fh2TracksLeadingPhiEta(0x0),
208 fh2TracksLeadingPhiPt(0x0),
209 fh2TracksLeadingJetPhiPt(0x0),
210 fh2JetsLeadingPhiPtW(0x0),
211 fh2TracksLeadingPhiPtW(0x0),
212 fh2TracksLeadingJetPhiPtW(0x0),
213 fh2NRecJetsPtRan(0x0),
215 fh2NConstLeadingPtRan(0x0),
220 fh2TracksLeadingJetPhiPtRan(0x0),
221 fh2TracksLeadingJetPhiPtWRan(0x0),
222 fh3CentvsRhoLeadingTrackPt(0x0),
223 fh3CentvsSigmaLeadingTrackPt(0x0),
224 fh3MultvsRhoLeadingTrackPt(0x0),
225 fh3MultvsSigmaLeadingTrackPt(0x0),
226 fh3CentvsRhoLeadingTrackPtQ1(0x0),
227 fh3CentvsRhoLeadingTrackPtQ2(0x0),
228 fh3CentvsRhoLeadingTrackPtQ3(0x0),
229 fh3CentvsRhoLeadingTrackPtQ4(0x0),
230 fh3CentvsSigmaLeadingTrackPtQ1(0x0),
231 fh3CentvsSigmaLeadingTrackPtQ2(0x0),
232 fh3CentvsSigmaLeadingTrackPtQ3(0x0),
233 fh3CentvsSigmaLeadingTrackPtQ4(0x0),
234 fh3MultvsRhoLeadingTrackPtQ1(0x0),
235 fh3MultvsRhoLeadingTrackPtQ2(0x0),
236 fh3MultvsRhoLeadingTrackPtQ3(0x0),
237 fh3MultvsRhoLeadingTrackPtQ4(0x0),
238 fh3MultvsSigmaLeadingTrackPtQ1(0x0),
239 fh3MultvsSigmaLeadingTrackPtQ2(0x0),
240 fh3MultvsSigmaLeadingTrackPtQ3(0x0),
241 fh3MultvsSigmaLeadingTrackPtQ4(0x0),
242 fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0),
243 fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0),
244 fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0),
245 fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0),
246 fh2PtGenPtSmeared(0x0),
248 fp1PtResolution(0x0),
255 for(int i = 0;i<3;i++){
256 fh1BiARandomCones[i] = 0;
257 fh1BiARandomConesRan[i] = 0;
259 for(int i = 0;i<kMaxCent;i++){
260 fh2JetsLeadingPhiPtC[i] = 0;
261 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
262 fh2TracksLeadingJetPhiPtC[i] = 0;
263 fh2TracksLeadingJetPhiPtWC[i] = 0;
267 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
268 AliAnalysisTaskSE(name),
272 fUseAODTrackInput(kFALSE),
273 fUseAODMCInput(kFALSE),
274 fUseBackgroundCalc(kFALSE),
275 fEventSelection(kFALSE),
276 fRequireVZEROAC(kFALSE),
277 fRequireTZEROvtx(kFALSE),
279 fFilterMaskBestPt(0),
282 fTrackTypeRec(kTrackUndef),
283 fTrackTypeGen(kTrackUndef),
285 fNSkipLeadingCone(0),
289 fTrackEtaWindow(0.9),
291 fApplySharedClusterCut(0),
294 fJetOutputMinPt(0.150),
295 fMaxTrackPtInJet(100.),
301 fStoreRhoLeadingTrackCorr(kFALSE),
303 fBackgroundBranch(""),
314 fUseTrPtResolutionSmearing(kFALSE),
315 fUseDiceEfficiency(kFALSE),
316 fDiceEfficiencyMinPt(-1.),
317 fUseTrPtResolutionFromOADB(kFALSE),
318 fUseTrEfficiencyFromOADB(kFALSE),
319 fPathTrPtResolution(""),
320 fPathTrEfficiency(""),
321 fChangeEfficiencyFraction(0.),
322 fEfficiencyFixed(1.),
324 fAlgorithm(fastjet::kt_algorithm),
325 fStrategy(fastjet::Best),
326 fRecombScheme(fastjet::BIpt_scheme),
327 fAreaType(fastjet::active_area),
329 fActiveAreaRepeats(1),
333 fTCARandomConesOut(0x0),
334 fTCARandomConesOutRan(0x0),
335 fAODJetBackgroundOut(0x0),
341 fh1PtHardTrials(0x0),
344 fh1NConstLeadingRec(0x0),
346 fh1PtJetsLeadingRecIn(0x0),
347 fh1PtJetConstRec(0x0),
348 fh1PtJetConstLeadingRec(0x0),
349 fh1PtTracksRecIn(0x0),
350 fh1PtTracksLeadingRecIn(0x0),
352 fh1NConstRecRan(0x0),
353 fh1PtJetsLeadingRecInRan(0x0),
354 fh1NConstLeadingRecRan(0x0),
355 fh1PtJetsRecInRan(0x0),
356 fh1PtTracksGenIn(0x0),
358 fh1CentralityPhySel(0x0),
360 fh1CentralitySelect(0x0),
365 fh2NRecTracksPt(0x0),
367 fh2NConstLeadingPt(0x0),
369 fh2LeadingJetPhiEta(0x0),
371 fh2LeadingJetEtaPt(0x0),
373 fh2LeadingTrackEtaPt(0x0),
374 fh2JetsLeadingPhiEta(0x0),
375 fh2JetsLeadingPhiPt(0x0),
376 fh2TracksLeadingPhiEta(0x0),
377 fh2TracksLeadingPhiPt(0x0),
378 fh2TracksLeadingJetPhiPt(0x0),
379 fh2JetsLeadingPhiPtW(0x0),
380 fh2TracksLeadingPhiPtW(0x0),
381 fh2TracksLeadingJetPhiPtW(0x0),
382 fh2NRecJetsPtRan(0x0),
384 fh2NConstLeadingPtRan(0x0),
389 fh2TracksLeadingJetPhiPtRan(0x0),
390 fh2TracksLeadingJetPhiPtWRan(0x0),
391 fh3CentvsRhoLeadingTrackPt(0x0),
392 fh3CentvsSigmaLeadingTrackPt(0x0),
393 fh3MultvsRhoLeadingTrackPt(0x0),
394 fh3MultvsSigmaLeadingTrackPt(0x0),
395 fh3CentvsRhoLeadingTrackPtQ1(0x0),
396 fh3CentvsRhoLeadingTrackPtQ2(0x0),
397 fh3CentvsRhoLeadingTrackPtQ3(0x0),
398 fh3CentvsRhoLeadingTrackPtQ4(0x0),
399 fh3CentvsSigmaLeadingTrackPtQ1(0x0),
400 fh3CentvsSigmaLeadingTrackPtQ2(0x0),
401 fh3CentvsSigmaLeadingTrackPtQ3(0x0),
402 fh3CentvsSigmaLeadingTrackPtQ4(0x0),
403 fh3MultvsRhoLeadingTrackPtQ1(0x0),
404 fh3MultvsRhoLeadingTrackPtQ2(0x0),
405 fh3MultvsRhoLeadingTrackPtQ3(0x0),
406 fh3MultvsRhoLeadingTrackPtQ4(0x0),
407 fh3MultvsSigmaLeadingTrackPtQ1(0x0),
408 fh3MultvsSigmaLeadingTrackPtQ2(0x0),
409 fh3MultvsSigmaLeadingTrackPtQ3(0x0),
410 fh3MultvsSigmaLeadingTrackPtQ4(0x0),
411 fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0),
412 fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0),
413 fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0),
414 fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0),
415 fh2PtGenPtSmeared(0x0),
417 fp1PtResolution(0x0),
424 for(int i = 0;i<3;i++){
425 fh1BiARandomCones[i] = 0;
426 fh1BiARandomConesRan[i] = 0;
428 for(int i = 0;i<kMaxCent;i++){
429 fh2JetsLeadingPhiPtC[i] = 0;
430 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
431 fh2TracksLeadingJetPhiPtC[i] = 0;
432 fh2TracksLeadingJetPhiPtWC[i] = 0;
434 DefineOutput(1, TList::Class());
439 Bool_t AliAnalysisTaskJetCluster::Notify()
442 // Implemented Notify() to read the cross sections
443 // and number of trials from pyxsec.root
448 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
452 // Create the output container
455 fRandom = new TRandom3(0);
461 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
465 if(fNonStdBranch.Length()!=0)
467 // only create the output branch if we have a name
468 // Create a new branch for jets...
469 // -> cleared in the UserExec....
470 // here we can also have the case that the brnaches are written to a separate file
473 fTCAJetsOut = new TClonesArray("AliAODJet", 0);
474 fTCAJetsOut->SetName(fNonStdBranch.Data());
475 AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
478 if(fJetTypes&kJetRan){
479 fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
480 fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
481 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
482 if( fEfficiencyFixed < 1.)
483 fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dEffFixed%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.)));
485 fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dFr%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
487 AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
490 if(fUseBackgroundCalc){
491 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
492 fAODJetBackgroundOut = new AliAODJetEventBackground();
493 fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
494 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
495 if( fEfficiencyFixed < 1.)
496 fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dEffFixed%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.)));
498 fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dFr%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
500 AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
503 // create the branch for the random cones with the same R
504 TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
505 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
506 if( fEfficiencyFixed < 1.)
507 cName = Form("%sDetector%d%dEffFixed%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.),fNSkipLeadingCone);
509 cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone);
513 if(!AODEvent()->FindListObject(cName.Data())){
514 fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
515 fTCARandomConesOut->SetName(cName.Data());
516 AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
519 // create the branch with the random for the random cones on the random event
520 if(fJetTypes&kRCRan){
521 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
522 if(!AODEvent()->FindListObject(cName.Data())){
523 fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
524 fTCARandomConesOutRan->SetName(cName.Data());
525 AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
530 if(fNonStdFile.Length()!=0){
532 // case that we have an AOD extension we need to fetch the jets from the extended output
533 // we identify the extension aod event by looking for the branchname
534 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
535 // case that we have an AOD extension we need can fetch the background maybe from the extended output
536 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
541 if(!fHistList)fHistList = new TList();
542 fHistList->SetOwner();
543 PostData(1, fHistList); // post data in any case once
545 Bool_t oldStatus = TH1::AddDirectoryStatus();
546 TH1::AddDirectory(kFALSE);
551 const Int_t nBinPt = 100;
552 Double_t binLimitsPt[nBinPt+1];
553 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
555 binLimitsPt[iPt] = 0.0;
558 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
562 const Int_t nBinPhi = 90;
563 Double_t binLimitsPhi[nBinPhi+1];
564 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
566 binLimitsPhi[iPhi] = -1.*TMath::Pi();
569 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
575 const Int_t nBinEta = 40;
576 Double_t binLimitsEta[nBinEta+1];
577 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
579 binLimitsEta[iEta] = -2.0;
582 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
586 const int nChMax = 5000;
588 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
589 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
591 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
592 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
595 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
596 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
598 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
599 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
600 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
601 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
604 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
605 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
606 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
608 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
609 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
610 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
611 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
612 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
613 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
614 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
615 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
616 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
617 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
619 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
620 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
621 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
623 fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
624 fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
625 fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
627 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
628 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
629 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
633 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
634 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
635 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
636 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
638 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
639 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);
640 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);
641 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);
645 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
646 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
647 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
648 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
650 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
651 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
652 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
653 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
655 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
656 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
657 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
658 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
662 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
663 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
664 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
665 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
666 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
667 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
668 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
669 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
670 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
671 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
672 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
673 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
675 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
676 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
677 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
678 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
680 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
681 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
682 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
683 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
685 if(fStoreRhoLeadingTrackCorr) {
686 fh3CentvsRhoLeadingTrackPt = new TH3F("fh3CentvsRhoLeadingTrackPt","centrality vs background density full event; centrality; #rho", 50,0.,100.,500,0.,250.,100,0.,100.);
687 fh3CentvsSigmaLeadingTrackPt = new TH3F("fh3CentvsSigmaLeadingTrackPt","centrality vs sigma full event; centrality; #sigma(#rho)", 50,0.,100.,50,0.,50.,100,0.,100.);
688 fh3MultvsRhoLeadingTrackPt = new TH3F("fh3MultvsRhoLeadingTrackPt","multiplicity vs background density full event; multiplicity; #rho", 100,0.,5000.,500,0.,250.,100,0.,100.);
689 fh3MultvsSigmaLeadingTrackPt = new TH3F("fh3MultvsSigmaLeadingTrackPt","multiplicity vs sigma full event; multiplicity; #sigma(#rho)", 100,0.,5000.,50,0.,50.,100,0.,100.);
692 fh3CentvsRhoLeadingTrackPtQ1 = new TH3F("fh3CentvsRhoLeadingTrackPtQ1","centrality vs background density Q1; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
693 fh3CentvsRhoLeadingTrackPtQ2 = new TH3F("fh3CentvsRhoLeadingTrackPtQ2","centrality vs background density Q2; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
694 fh3CentvsRhoLeadingTrackPtQ3 = new TH3F("fh3CentvsRhoLeadingTrackPtQ3","centrality vs background density Q3; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
695 fh3CentvsRhoLeadingTrackPtQ4 = new TH3F("fh3CentvsRhoLeadingTrackPtQ4","centrality vs background density Q4; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
697 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.);
698 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.);
699 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.);
700 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.);
702 fh3MultvsRhoLeadingTrackPtQ1 = new TH3F("fh3MultvsRhoLeadingTrackPtQ1","multiplicity vs background density Q1; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
703 fh3MultvsRhoLeadingTrackPtQ2 = new TH3F("fh3MultvsRhoLeadingTrackPtQ2","multiplicity vs background density Q2; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
704 fh3MultvsRhoLeadingTrackPtQ3 = new TH3F("fh3MultvsRhoLeadingTrackPtQ3","multiplicity vs background density Q3; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
705 fh3MultvsRhoLeadingTrackPtQ4 = new TH3F("fh3MultvsRhoLeadingTrackPtQ4","multiplicity vs background density Q4; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
707 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.);
708 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.);
709 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.);
710 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.);
713 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.);
714 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.);
715 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.);
716 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.);
718 fHistList->Add(fh3CentvsRhoLeadingTrackPt);
719 fHistList->Add(fh3CentvsSigmaLeadingTrackPt);
720 fHistList->Add(fh3MultvsRhoLeadingTrackPt);
721 fHistList->Add(fh3MultvsSigmaLeadingTrackPt);
723 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ1);
724 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ2);
725 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ3);
726 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ4);
728 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ1);
729 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ2);
730 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ3);
731 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ4);
733 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ1);
734 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ2);
735 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ3);
736 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ4);
738 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ1);
739 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ2);
740 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ3);
741 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ4);
743 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ1);
744 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ2);
745 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ3);
746 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ4);
750 //Detector level effects histos
751 fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
753 fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
754 fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
756 fHistList->Add(fh2PtGenPtSmeared);
757 fHistList->Add(fp1Efficiency);
758 fHistList->Add(fp1PtResolution);
760 if(fNRandomCones>0&&fUseBackgroundCalc){
761 for(int i = 0;i<3;i++){
762 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
763 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
767 for(int i = 0;i < kMaxCent;i++){
768 fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
769 fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
770 fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
771 fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
774 const Int_t saveLevel = 3; // large save level more histos
776 fHistList->Add(fh1Xsec);
777 fHistList->Add(fh1Trials);
779 fHistList->Add(fh1NJetsRec);
780 fHistList->Add(fh1NConstRec);
781 fHistList->Add(fh1NConstLeadingRec);
782 fHistList->Add(fh1PtJetsRecIn);
783 fHistList->Add(fh1PtJetsLeadingRecIn);
784 fHistList->Add(fh1PtTracksRecIn);
785 fHistList->Add(fh1PtTracksLeadingRecIn);
786 fHistList->Add(fh1PtJetConstRec);
787 fHistList->Add(fh1PtJetConstLeadingRec);
788 fHistList->Add(fh1NJetsRecRan);
789 fHistList->Add(fh1NConstRecRan);
790 fHistList->Add(fh1PtJetsLeadingRecInRan);
791 fHistList->Add(fh1NConstLeadingRecRan);
792 fHistList->Add(fh1PtJetsRecInRan);
793 fHistList->Add(fh1Nch);
794 fHistList->Add(fh1Centrality);
795 fHistList->Add(fh1CentralitySelect);
796 fHistList->Add(fh1CentralityPhySel);
797 fHistList->Add(fh1Z);
798 fHistList->Add(fh1ZSelect);
799 fHistList->Add(fh1ZPhySel);
800 if(fNRandomCones>0&&fUseBackgroundCalc){
801 for(int i = 0;i<3;i++){
802 fHistList->Add(fh1BiARandomCones[i]);
803 fHistList->Add(fh1BiARandomConesRan[i]);
806 for(int i = 0;i < kMaxCent;i++){
807 fHistList->Add(fh2JetsLeadingPhiPtC[i]);
808 fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
809 fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
810 fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
813 fHistList->Add(fh2NRecJetsPt);
814 fHistList->Add(fh2NRecTracksPt);
815 fHistList->Add(fh2NConstPt);
816 fHistList->Add(fh2NConstLeadingPt);
817 fHistList->Add(fh2PtNch);
818 fHistList->Add(fh2PtNchRan);
819 fHistList->Add(fh2PtNchN);
820 fHistList->Add(fh2PtNchNRan);
821 fHistList->Add(fh2JetPhiEta);
822 fHistList->Add(fh2LeadingJetPhiEta);
823 fHistList->Add(fh2JetEtaPt);
824 fHistList->Add(fh2LeadingJetEtaPt);
825 fHistList->Add(fh2TrackEtaPt);
826 fHistList->Add(fh2LeadingTrackEtaPt);
827 fHistList->Add(fh2JetsLeadingPhiEta );
828 fHistList->Add(fh2JetsLeadingPhiPt);
829 fHistList->Add(fh2TracksLeadingPhiEta);
830 fHistList->Add(fh2TracksLeadingPhiPt);
831 fHistList->Add(fh2TracksLeadingJetPhiPt);
832 fHistList->Add(fh2JetsLeadingPhiPtW);
833 fHistList->Add(fh2TracksLeadingPhiPtW);
834 fHistList->Add(fh2TracksLeadingJetPhiPtW);
835 fHistList->Add(fh2NRecJetsPtRan);
836 fHistList->Add(fh2NConstPtRan);
837 fHistList->Add(fh2NConstLeadingPtRan);
838 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
839 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
842 // =========== Switch on Sumw2 for all histos ===========
843 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
844 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
849 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
852 TH1::AddDirectory(oldStatus);
855 void AliAnalysisTaskJetCluster::LocalInit()
861 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
863 if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB();
864 if(fUseTrEfficiencyFromOADB) LoadTrEfficiencyRootFileFromOADB();
867 FitMomentumResolution();
871 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
874 // handle and reset the output jet branch
876 if(fTCAJetsOut)fTCAJetsOut->Delete();
877 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
878 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
879 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
880 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
882 AliAODJetEventBackground* externalBackground = 0;
883 if(!externalBackground&&fBackgroundBranch.Length()){
884 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
885 if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
886 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
889 // Execute analysis for current event
891 AliESDEvent *fESD = 0;
892 if(fUseAODTrackInput){
893 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
895 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
901 // assume that the AOD is in the general output...
904 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
908 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
912 //Check if information is provided detector level effects
913 if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE;
914 if( fEfficiencyFixed < 1. ) {
915 if (!fUseDiceEfficiency)
916 fUseDiceEfficiency = 1; // 1 is the default; 2 can be set by user
919 if(!fhEffH1 || !fhEffH2 || !fhEffH3 ) fUseDiceEfficiency = kFALSE;
922 Bool_t selectEvent = false;
923 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
929 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
930 TString vtxTitle(vtxAOD->GetTitle());
931 zVtx = vtxAOD->GetZ();
933 cent = fAOD->GetHeader()->GetCentrality();
934 if(cent<10)cenClass = 0;
935 else if(cent<30)cenClass = 1;
936 else if(cent<50)cenClass = 2;
937 else if(cent<80)cenClass = 3;
938 if(physicsSelection){
939 fh1CentralityPhySel->Fill(cent);
940 fh1ZPhySel->Fill(zVtx);
944 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
945 Float_t yvtx = vtxAOD->GetY();
946 Float_t xvtx = vtxAOD->GetX();
947 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
948 if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
949 if(physicsSelection){
955 if(cent<fCentCutLo||cent>fCentCutUp){
967 const AliAODVZERO *vzero = fAOD->GetVZEROData();
969 if((vzero->GetTriggerChargeA()>0)&&(vzero->GetTriggerChargeC()>0)){
974 const AliAODTZERO *tzero = fAOD->GetTZEROData();
976 if(TMath::Abs(tzero->GetT0VertexRaw())<100){
981 if(fRequireVZEROAC&&fRequireTZEROvtx)selectEvent = selectEvent&&V0&&T0;
982 else if(fRequireTZEROvtx)selectEvent = selectEvent&&T0;
983 else if(fRequireVZEROAC)selectEvent = selectEvent&&V0;
987 PostData(1, fHistList);
990 fh1Centrality->Fill(cent);
992 fh1Trials->Fill("#sum{ntrials}",1);
995 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
997 // ==== General variables needed
1001 // we simply fetch the tracks/mc particles as a list of AliVParticles
1006 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
1007 Float_t nCh = recParticles.GetEntries();
1008 Float_t nGen=genParticles.GetEntries();
1010 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
1011 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
1012 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
1014 // find the jets....
1016 vector<fastjet::PseudoJet> inputParticlesRec;
1017 vector<fastjet::PseudoJet> inputParticlesRecRan;
1019 // Generate the random cones
1021 AliAODJet vTmpRan(1,0,0,1);
1022 for(int i = 0; i < recParticles.GetEntries(); i++){
1023 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1025 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1026 // we take total momentum here
1028 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
1029 //Add particles to fastjet in case we are not running toy model
1030 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
1031 jInp.set_user_index(i);
1032 inputParticlesRec.push_back(jInp);
1034 else if(fUseDiceEfficiency) {
1036 // Dice to decide if particle is kept or not - toy model for efficiency
1038 Double_t sumEff = 0.;
1040 Double_t eff[3] = {0.};
1043 Double_t rnd = fRandom->Uniform(1.);
1044 if( fEfficiencyFixed<1. ) {
1045 sumEff = fEfficiencyFixed;
1046 if (fUseDiceEfficiency == 2) {
1047 sumEff = (nCh - fEfficiencyFixed*nGen) / nCh; // rescale eff; fEfficiencyFixed is wrt to nGen, but dicing is fraction of nCh
1052 Double_t pTtmp = pT;
1053 if(pT>10.) pTtmp = 10.;
1054 Double_t eff1 = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
1055 Double_t eff2 = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
1056 Double_t eff3 = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));
1058 //Sort efficiencies from large to small
1059 if(eff1>eff2 && eff1>eff3) {
1074 else if(eff2>eff1 && eff2>eff3) {
1089 else if(eff3>eff1 && eff3>eff2) {
1105 sumEff = eff[0]+eff[1]+eff[2];
1107 fp1Efficiency->Fill(vp->Pt(),sumEff);
1108 if(rnd>sumEff && pT > fDiceEfficiencyMinPt) continue;
1110 if(fUseTrPtResolutionSmearing) {
1111 //Smear momentum of generated particle
1112 Double_t smear = 1.;
1113 //Select hybrid track category
1115 smear = GetMomentumSmearing(cat[2],pT);
1116 else if(rnd<=(eff[2]+eff[1]))
1117 smear = GetMomentumSmearing(cat[1],pT);
1119 smear = GetMomentumSmearing(cat[0],pT);
1121 fp1PtResolution->Fill(vp->Pt(),smear);
1123 Double_t sigma = vp->Pt()*smear;
1124 Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
1126 Double_t phi = vp->Phi();
1127 Double_t theta = 2.*TMath::ATan(TMath::Exp(-1.*(vp->Eta())));
1128 Double_t pX = pTrec * TMath::Cos(phi);
1129 Double_t pY = pTrec * TMath::Sin(phi);
1130 Double_t pZ = pTrec/TMath::Tan(theta);
1131 Double_t p=TMath::Sqrt(pTrec*pTrec+pZ*pZ);
1133 fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);
1135 fastjet::PseudoJet jInp(pX,pY,pZ,p);
1136 jInp.set_user_index(i);
1137 inputParticlesRec.push_back(jInp);
1141 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
1142 jInp.set_user_index(i);
1143 inputParticlesRec.push_back(jInp);
1149 // the randomized input changes eta and phi, but keeps the p_T
1150 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1151 Double_t pT = vp->Pt();
1152 Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
1153 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
1155 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
1156 Double_t pZ = pT/TMath::Tan(theta);
1158 Double_t pX = pT * TMath::Cos(phi);
1159 Double_t pY = pT * TMath::Sin(phi);
1160 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
1161 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
1163 jInpRan.set_user_index(i);
1164 inputParticlesRecRan.push_back(jInpRan);
1165 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
1168 // fill the tref array, only needed when we write out jets
1171 fRef->Delete(); // make sure to delete before placement new...
1172 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
1173 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
1176 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ...
1180 if(inputParticlesRec.size()==0){
1181 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
1182 PostData(1, fHistList);
1187 // employ setters for these...
1190 // now create the object that holds info about ghosts
1192 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
1193 // reduce CPU time...
1195 fActiveAreaRepeats = 0;
1199 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
1200 fastjet::AreaType areaType = fastjet::active_area;
1201 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
1202 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
1203 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
1205 //range where to compute background
1206 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
1208 phiMax = 2*TMath::Pi();
1209 rapMax = fGhostEtamax - fRparam;
1210 rapMin = - fGhostEtamax + fRparam;
1211 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
1214 const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
1215 const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
1218 fh1NJetsRec->Fill(sortedJets.size());
1220 // loop over all jets an fill information, first on is the leading jet
1222 Int_t nRecOver = inclusiveJets.size();
1223 Int_t nRec = inclusiveJets.size();
1224 if(inclusiveJets.size()>0){
1225 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
1226 Double_t area = clustSeq.area(sortedJets[0]);
1227 leadingJet.SetEffArea(area,0);
1228 Float_t pt = leadingJet.Pt();
1229 Int_t nAodOutJets = 0;
1230 Int_t nAodOutTracks = 0;
1231 AliAODJet *aodOutJet = 0;
1234 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
1235 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
1236 while(pt<ptCut&&iCount<nRec){
1240 pt = sortedJets[iCount].perp();
1243 if(nRecOver<=0)break;
1244 fh2NRecJetsPt->Fill(ptCut,nRecOver);
1246 Float_t phi = leadingJet.Phi();
1247 if(phi<0)phi+=TMath::Pi()*2.;
1248 Float_t eta = leadingJet.Eta();
1250 if(externalBackground){
1251 // carefull has to be filled in a task before
1252 // todo, ReArrange to the botom
1253 pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
1255 pt = leadingJet.Pt() - pTback;
1256 // correlation of leading jet with tracks
1257 TIterator *recIter = recParticles.MakeIterator();
1259 AliVParticle *tmpRecTrack = 0;
1260 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
1261 Float_t tmpPt = tmpRecTrack->Pt();
1263 Float_t tmpPhi = tmpRecTrack->Phi();
1264 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1265 Float_t dPhi = phi - tmpPhi;
1266 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1267 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1268 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
1269 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
1271 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
1272 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1278 TLorentzVector vecareab;
1279 for(int j = 0; j < nRec;j++){
1280 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
1283 Float_t tmpPt = tmpRec.Pt();
1285 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
1286 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
1287 aodOutJet->GetRefTracks()->Clear();
1288 Double_t area1 = clustSeq.area(sortedJets[j]);
1289 aodOutJet->SetEffArea(area1,0);
1290 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
1291 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
1292 aodOutJet->SetVectorAreaCharged(&vecareab);
1296 Float_t tmpPtBack = 0;
1297 if(externalBackground){
1298 // carefull has to be filled in a task before
1299 // todo, ReArrange to the botom
1300 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
1302 tmpPt = tmpPt - tmpPtBack;
1303 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
1305 fh1PtJetsRecIn->Fill(tmpPt);
1306 // Fill Spectra with constituentsemacs
1307 const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
1309 fh1NConstRec->Fill(constituents.size());
1310 fh2PtNch->Fill(nCh,tmpPt);
1311 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
1312 fh2NConstPt->Fill(tmpPt,constituents.size());
1313 // loop over constiutents and fill spectrum
1315 AliVParticle *partLead = 0;
1316 Float_t ptLead = -1;
1318 for(unsigned int ic = 0; ic < constituents.size();ic++){
1319 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1321 fh1PtJetConstRec->Fill(part->Pt());
1323 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1324 if(part->Pt()>fMaxTrackPtInJet){
1325 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1328 if(part->Pt()>ptLead){
1329 ptLead = part->Pt();
1332 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1335 AliAODTrack *aodT = 0;
1338 //set pT of leading constituent of jet
1339 aodOutJet->SetPtLeading(partLead->Pt());
1340 aodT = dynamic_cast<AliAODTrack*>(partLead);
1342 if(aodT->TestFilterBit(fFilterMaskBestPt)){
1343 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
1350 Float_t tmpPhi = tmpRec.Phi();
1351 Float_t tmpEta = tmpRec.Eta();
1352 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1354 fh1PtJetsLeadingRecIn->Fill(tmpPt);
1355 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1356 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1357 fh1NConstLeadingRec->Fill(constituents.size());
1358 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1361 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1362 fh2JetEtaPt->Fill(tmpEta,tmpPt);
1363 Float_t dPhi = phi - tmpPhi;
1364 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1365 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1366 Float_t dEta = eta - tmpRec.Eta();
1367 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1368 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1370 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1371 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1373 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1374 }// loop over reconstructed jets
1379 // Add the random cones...
1380 if(fNRandomCones>0&&fTCARandomConesOut){
1381 // create a random jet within the acceptance
1382 Double_t etaMax = fTrackEtaWindow - fRparam;
1385 Double_t pTC = 1; // small number
1386 for(int ir = 0;ir < fNRandomCones;ir++){
1387 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1388 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1390 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1391 Double_t pZC = pTC/TMath::Tan(thetaC);
1392 Double_t pXC = pTC * TMath::Cos(phiC);
1393 Double_t pYC = pTC * TMath::Sin(phiC);
1394 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1395 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
1397 for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1398 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1399 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1404 // test for overlap with previous cones to avoid double counting
1405 for(int iic = 0;iic<ir;iic++){
1406 AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1408 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1415 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1416 tmpRecC.SetPtLeading(-1.);
1417 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1418 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1419 }// loop over random cones creation
1422 // loop over the reconstructed particles and add up the pT in the random cones
1423 // maybe better to loop over randomized particles not in the real jets...
1424 // but this by definition brings dow average energy in the whole event
1425 AliAODJet vTmpRanR(1,0,0,1);
1426 for(int i = 0; i < recParticles.GetEntries(); i++){
1427 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1428 if(fTCARandomConesOut){
1429 for(int ir = 0;ir < fNRandomCones;ir++){
1430 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1431 if(jC&&jC->DeltaR(vp)<fRparam){
1432 if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1433 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1434 if(vp->Pt() > jC->GetPtLeading()) jC->SetPtLeading(vp->Pt());
1437 }// add up energy in cone
1439 // the randomized input changes eta and phi, but keeps the p_T
1440 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1441 Double_t pTR = vp->Pt();
1442 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1443 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1445 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1446 Double_t pZR = pTR/TMath::Tan(thetaR);
1448 Double_t pXR = pTR * TMath::Cos(phiR);
1449 Double_t pYR = pTR * TMath::Sin(phiR);
1450 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1451 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1452 if(fTCARandomConesOutRan){
1453 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1454 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1455 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1456 if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1457 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1458 if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt());
1463 }// loop over recparticles
1465 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1466 if(fTCARandomConesOut){
1467 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1468 // rescale the momentum vectors for the random cones
1470 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1472 Double_t etaC = rC->Eta();
1473 Double_t phiC = rC->Phi();
1474 // massless jet, unit vector
1475 pTC = rC->ChargedBgEnergy();
1476 if(pTC<=0)pTC = 0.001; // for almost empty events
1477 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1478 Double_t pZC = pTC/TMath::Tan(thetaC);
1479 Double_t pXC = pTC * TMath::Cos(phiC);
1480 Double_t pYC = pTC * TMath::Sin(phiC);
1481 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1482 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1483 rC->SetBgEnergy(0,0);
1484 rC->SetEffArea(jetArea,0);
1488 if(fTCARandomConesOutRan){
1489 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1490 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1493 Double_t etaC = rC->Eta();
1494 Double_t phiC = rC->Phi();
1495 // massless jet, unit vector
1496 pTC = rC->ChargedBgEnergy();
1497 if(pTC<=0)pTC = 0.001;// for almost empty events
1498 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1499 Double_t pZC = pTC/TMath::Tan(thetaC);
1500 Double_t pXC = pTC * TMath::Cos(phiC);
1501 Double_t pYC = pTC * TMath::Sin(phiC);
1502 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1503 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1504 rC->SetBgEnergy(0,0);
1505 rC->SetEffArea(jetArea,0);
1509 }// if(fNRandomCones
1511 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1513 if(fAODJetBackgroundOut){
1514 vector<fastjet::PseudoJet> jets2=sortedJets;
1515 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1519 Double_t meanarea1=0.;
1522 Double_t meanarea2=0.;
1526 Double_t meanarea4=0.;
1528 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1529 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1531 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1532 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1533 clustSeq.get_median_rho_and_sigma(sortedJets, range, true, bkg4, sigma4, meanarea4, true);
1534 fAODJetBackgroundOut->SetBackground(3,bkg4,sigma4,meanarea4);
1536 //////////////////////
1538 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1539 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1540 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1541 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1546 if(fStoreRhoLeadingTrackCorr) {
1547 vector<fastjet::PseudoJet> jets3=sortedJets;
1548 if(jets3.size()>2) jets3.erase(jets3.begin(),jets3.begin()+2);
1552 Double_t meanarea2=0.;
1554 clustSeq.get_median_rho_and_sigma(jets3, range, false, bkg2, sigma2, meanarea2, true);
1555 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1557 //Get leading track in event
1558 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1560 fh3CentvsRhoLeadingTrackPt->Fill(cent,bkg2,leading->Pt());
1561 fh3CentvsSigmaLeadingTrackPt->Fill(cent,sigma2,leading->Pt());
1562 fh3MultvsRhoLeadingTrackPt->Fill(nCh,bkg2,leading->Pt());
1563 fh3MultvsSigmaLeadingTrackPt->Fill(nCh,sigma2,leading->Pt());
1566 //Calculate rho with 4-vector area method for different quadrants with respect to the leading track in the event
1567 //exclude 2 leading jets from event
1568 //Quadrant 1: |phi_leadingTrack - phi_bkgCluster| < pi/2./2. - R (Near side to leading track)
1569 //Quadrant 2: pi/2 - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < pi/2 + (pi/2./2. - R) (Orthogonal to leading track)
1570 //Quadrant 3: pi - (pi/2/2 - R) < |phi_leadingTrack - phi_bkgCluster| < pi + (pi/2./2. - R) (Away side to leading track)
1571 //Quadrant 4: 3/2*pi - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < 3/2*pi + (pi/2./2. - R) (Orthogonal to leading track)
1573 //Assuming R=0.2 for background clusters
1575 Double_t bkg2Q[4] = {0.};
1576 Double_t sigma2Q[4] = {0.};
1577 Double_t meanarea2Q[4] = {0.};
1579 //Get phi of leading track in event
1580 Float_t phiLeadingTrack = leading->Phi();
1581 Float_t phiStep = (TMath::Pi()/2./2. - 0.2);
1583 //Quadrant1 - near side
1584 phiMin = phiLeadingTrack - phiStep;
1585 phiMax = phiLeadingTrack + phiStep;
1586 fastjet::RangeDefinition rangeQ1(rapMin,rapMax, phiMin, phiMax);
1587 clustSeq.get_median_rho_and_sigma(jets3, rangeQ1, false, bkg2Q[0], sigma2Q[0], meanarea2Q[0], true);
1589 //Quadrant2 - orthogonal
1590 Double_t phiQ2 = phiLeadingTrack + TMath::Pi()/2.;
1591 if(phiQ2>TMath::TwoPi()) phiQ2 = phiQ2 - TMath::TwoPi();
1592 phiMin = phiQ2 - phiStep;
1593 phiMax = phiQ2 + phiStep;
1594 fastjet::RangeDefinition rangeQ2(rapMin,rapMax, phiMin, phiMax);
1595 clustSeq.get_median_rho_and_sigma(jets3, rangeQ2, false, bkg2Q[1], sigma2Q[1], meanarea2Q[1], true);
1597 //Quadrant3 - away side
1598 Double_t phiQ3 = phiLeadingTrack + TMath::Pi();
1599 if(phiQ3>TMath::TwoPi()) phiQ3 = phiQ3 - TMath::TwoPi();
1600 phiMin = phiQ3 - phiStep;
1601 phiMax = phiQ3 + phiStep;
1602 fastjet::RangeDefinition rangeQ3(rapMin,rapMax, phiMin, phiMax);
1603 clustSeq.get_median_rho_and_sigma(jets3, rangeQ3, false, bkg2Q[2], sigma2Q[2], meanarea2Q[2], true);
1605 //Quadrant4 - othogonal
1606 Double_t phiQ4 = phiLeadingTrack + 3./2.*TMath::Pi();
1607 if(phiQ4>TMath::TwoPi()) phiQ4 = phiQ4 - TMath::TwoPi();
1608 phiMin = phiQ4 - phiStep;
1609 phiMax = phiQ4 + phiStep;
1610 fastjet::RangeDefinition rangeQ4(rapMin,rapMax, phiMin, phiMax);
1611 clustSeq.get_median_rho_and_sigma(jets3, rangeQ4, false, bkg2Q[3], sigma2Q[3], meanarea2Q[3], true);
1613 fh3CentvsRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0],leading->Pt());
1614 fh3CentvsRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1],leading->Pt());
1615 fh3CentvsRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2],leading->Pt());
1616 fh3CentvsRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3],leading->Pt());
1618 fh3CentvsSigmaLeadingTrackPtQ1->Fill(cent,sigma2Q[0],leading->Pt());
1619 fh3CentvsSigmaLeadingTrackPtQ2->Fill(cent,sigma2Q[1],leading->Pt());
1620 fh3CentvsSigmaLeadingTrackPtQ3->Fill(cent,sigma2Q[2],leading->Pt());
1621 fh3CentvsSigmaLeadingTrackPtQ4->Fill(cent,sigma2Q[3],leading->Pt());
1623 fh3MultvsRhoLeadingTrackPtQ1->Fill(nCh,bkg2Q[0],leading->Pt());
1624 fh3MultvsRhoLeadingTrackPtQ2->Fill(nCh,bkg2Q[1],leading->Pt());
1625 fh3MultvsRhoLeadingTrackPtQ3->Fill(nCh,bkg2Q[2],leading->Pt());
1626 fh3MultvsRhoLeadingTrackPtQ4->Fill(nCh,bkg2Q[3],leading->Pt());
1628 fh3MultvsSigmaLeadingTrackPtQ1->Fill(nCh,sigma2Q[0],leading->Pt());
1629 fh3MultvsSigmaLeadingTrackPtQ2->Fill(nCh,sigma2Q[1],leading->Pt());
1630 fh3MultvsSigmaLeadingTrackPtQ3->Fill(nCh,sigma2Q[2],leading->Pt());
1631 fh3MultvsSigmaLeadingTrackPtQ4->Fill(nCh,sigma2Q[3],leading->Pt());
1633 fh3CentvsDeltaRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0]-bkg2,leading->Pt());
1634 fh3CentvsDeltaRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1]-bkg2,leading->Pt());
1635 fh3CentvsDeltaRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2]-bkg2,leading->Pt());
1636 fh3CentvsDeltaRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3]-bkg2,leading->Pt());
1644 // fill track information
1645 Int_t nTrackOver = recParticles.GetSize();
1646 // do the same for tracks and jets
1649 TIterator *recIter = recParticles.MakeIterator();
1650 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1651 Float_t pt = tmpRec->Pt();
1653 // Printf("Leading track p_t %3.3E",pt);
1654 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1655 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1656 while(pt<ptCut&&tmpRec){
1658 tmpRec = (AliVParticle*)(recIter->Next());
1663 if(nTrackOver<=0)break;
1664 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1668 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1669 Float_t phi = leading->Phi();
1670 if(phi<0)phi+=TMath::Pi()*2.;
1671 Float_t eta = leading->Eta();
1673 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1674 Float_t tmpPt = tmpRec->Pt();
1675 Float_t tmpEta = tmpRec->Eta();
1676 fh1PtTracksRecIn->Fill(tmpPt);
1677 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1678 if(tmpRec==leading){
1679 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1680 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1684 Float_t tmpPhi = tmpRec->Phi();
1686 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1687 Float_t dPhi = phi - tmpPhi;
1688 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1689 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1690 Float_t dEta = eta - tmpRec->Eta();
1691 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1692 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1693 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1698 // find the random jets
1700 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1702 // fill the jet information from random track
1703 const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1704 const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1706 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1708 // loop over all jets an fill information, first on is the leading jet
1710 Int_t nRecOverRan = inclusiveJetsRan.size();
1711 Int_t nRecRan = inclusiveJetsRan.size();
1713 if(inclusiveJetsRan.size()>0){
1714 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1715 Float_t pt = leadingJet.Pt();
1718 TLorentzVector vecarearanb;
1720 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1721 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1722 while(pt<ptCut&&iCount<nRecRan){
1726 pt = sortedJetsRan[iCount].perp();
1729 if(nRecOverRan<=0)break;
1730 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1732 Float_t phi = leadingJet.Phi();
1733 if(phi<0)phi+=TMath::Pi()*2.;
1734 pt = leadingJet.Pt();
1736 // correlation of leading jet with random tracks
1738 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1740 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1742 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1743 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1744 Float_t dPhi = phi - tmpPhi;
1745 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1746 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1747 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1748 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1751 Int_t nAodOutJetsRan = 0;
1752 AliAODJet *aodOutJetRan = 0;
1753 for(int j = 0; j < nRecRan;j++){
1754 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1755 Float_t tmpPt = tmpRec.Pt();
1756 fh1PtJetsRecInRan->Fill(tmpPt);
1757 // Fill Spectra with constituents
1758 const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1759 fh1NConstRecRan->Fill(constituents.size());
1760 fh2NConstPtRan->Fill(tmpPt,constituents.size());
1761 fh2PtNchRan->Fill(nCh,tmpPt);
1762 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1765 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1766 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1767 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1768 aodOutJetRan->GetRefTracks()->Clear();
1769 aodOutJetRan->SetEffArea(arearan,0);
1770 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1771 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1772 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1777 Float_t tmpPhi = tmpRec.Phi();
1778 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1781 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1782 fh1NConstLeadingRecRan->Fill(constituents.size());
1783 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1789 if(fAODJetBackgroundOut){
1792 Double_t meanarea3=0.;
1793 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1794 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1795 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1797 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1798 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1807 // do the event selection if activated
1808 if(fJetTriggerPtCut>0){
1809 bool select = false;
1810 Float_t minPt = fJetTriggerPtCut;
1812 // hard coded for now ...
1813 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1814 if(cent<10)minPt = 50;
1815 else if(cent<30)minPt = 42;
1816 else if(cent<50)minPt = 28;
1817 else if(cent<80)minPt = 18;
1820 if(externalBackground)rho = externalBackground->GetBackground(2);
1822 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1823 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1824 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1833 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1834 fh1CentralitySelect->Fill(cent);
1835 fh1ZSelect->Fill(zVtx);
1836 aodH->SetFillAOD(kTRUE);
1840 if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1841 if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1842 if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1843 if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1845 PostData(1, fHistList);
1848 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1851 // Terminate analysis
1853 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1855 if(fMomResH1Fit) delete fMomResH1Fit;
1856 if(fMomResH2Fit) delete fMomResH2Fit;
1857 if(fMomResH3Fit) delete fMomResH3Fit;
1862 Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1865 // get list of tracks/particles for different types
1868 if(fDebug>2) Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1871 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly || type==kTrackAODMCextra || type==kTrackAODMCextraonly){
1873 if(type!=kTrackAODextraonly && type!=kTrackAODMCextraonly) {
1874 AliAODEvent *aod = 0;
1875 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1876 else aod = AODEvent();
1878 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1882 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1883 AliAODTrack *tr = aod->GetTrack(it);
1884 Bool_t bGood = false;
1885 if(fFilterType == 0)bGood = true;
1886 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1887 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1888 if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1889 if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1892 if(fRequireITSRefit){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
1893 if (fApplySharedClusterCut) {
1894 Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
1895 if (frac > 0.4) continue;
1897 if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1898 if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1901 if(tr->Pt()<fTrackPtCut){
1902 if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1905 if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1910 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1911 AliAODEvent *aod = 0;
1912 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1913 else aod = AODEvent();
1918 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1919 if(!aodExtraTracks)return iCount;
1920 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1921 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1922 if (!track) continue;
1924 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1925 if(!trackAOD)continue;
1926 Bool_t bGood = false;
1927 if(fFilterType == 0)bGood = true;
1928 else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1929 else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1930 if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1931 if(fRequireITSRefit){if((trackAOD->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
1932 if (fApplySharedClusterCut) {
1933 Double_t frac = Double_t(trackAOD->GetTPCnclsS()) /Double_t(trackAOD->GetTPCncls());
1934 if (frac > 0.4) continue;
1938 if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1939 if(trackAOD->Pt()<fTrackPtCut) continue;
1940 if(fDebug) printf("pt extra track %.2f \n", trackAOD->Pt());
1941 list->Add(trackAOD);
1946 if(type==kTrackAODMCextra || type==kTrackAODMCextraonly) { //embed generator level particles
1947 AliAODEvent *aod = 0;
1948 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1949 else aod = AODEvent();
1953 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraMCparticles"));
1954 if(!aodExtraTracks)return iCount;
1955 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1956 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1957 AliAODMCParticle *partmc = dynamic_cast<AliAODMCParticle*> ((*aodExtraTracks)[it]);
1959 if(fDebug>10) printf("track %d does not exist\n",it);
1963 if(!partmc) continue;
1964 if(!partmc->IsPhysicalPrimary())continue;
1966 if (track->Pt()<fTrackPtCut) {
1967 if(fDebug>10) printf("track %d has too low pt %.2f\n",it,track->Pt());
1972 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*>((*aodExtraTracks)[it]);//(track);
1975 if(fDebug>10) printf("trackAOD %d does not exist\n",it);
1979 trackAOD->SetFlags(AliESDtrack::kEmbedded);
1980 trackAOD->SetFilterMap(fFilterMask);
1982 if(fDebug>10) printf("pt extra track %.2f \n", track->Pt());
1984 if(TMath::Abs(track->Eta())>fTrackEtaWindow) continue;
1985 if(track->Pt()<fTrackPtCut) continue;
1993 else if (type == kTrackKineAll||type == kTrackKineCharged){
1994 AliMCEvent* mcEvent = MCEvent();
1995 if(!mcEvent)return iCount;
1996 // we want to have alivpartilces so use get track
1997 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1998 if(!mcEvent->IsPhysicalPrimary(it))continue;
1999 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
2000 if(type == kTrackKineAll){
2001 if(part->Pt()<fTrackPtCut)continue;
2005 else if(type == kTrackKineCharged){
2006 if(part->Particle()->GetPDG()->Charge()==0)continue;
2007 if(part->Pt()<fTrackPtCut)continue;
2013 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
2014 AliAODEvent *aod = 0;
2015 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
2016 else aod = AODEvent();
2017 if(!aod)return iCount;
2018 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
2019 if(!tca)return iCount;
2020 for(int it = 0;it < tca->GetEntriesFast();++it){
2021 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
2022 if(!part->IsPhysicalPrimary())continue;
2023 if(type == kTrackAODMCAll){
2024 if(part->Pt()<fTrackPtCut)continue;
2028 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
2029 if(part->Charge()==0)continue;
2030 if(part->Pt()<fTrackPtCut)continue;
2031 if(kTrackAODMCCharged){
2035 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
2042 else if (type == kTrackAODMCHF){
2044 AliAODEvent *aod = 0;
2045 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
2046 else aod = AODEvent();
2047 if(!aod)return iCount;
2048 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
2049 if(!tca)return iCount;
2050 for(int it = 0;it < tca->GetEntriesFast();++it){
2051 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
2053 int partpdg= part->PdgCode();
2054 if(!part->IsPhysicalPrimary() && !IsBMeson(partpdg) && !IsDMeson(partpdg) )continue;
2056 if (IsBMeson(partpdg) || IsDMeson(partpdg)) {
2057 iCount+= AddDaughters( list , part,tca);
2061 if(part->Pt()<fTrackPtCut) continue;
2062 if(TMath::Abs(part->Eta())>fTrackEtaWindow) continue;
2063 if(part->Charge()==0) continue;
2065 if((part->Pt()>=fTrackPtCut) && (TMath::Abs(part->Eta())<=fTrackEtaWindow) && (part->Charge()!=0))list->Add(part);
2075 Int_t AliAnalysisTaskJetCluster::AddDaughters(TList * list, AliAODMCParticle *part, TClonesArray * tca){
2077 Int_t nDaugthers = part->GetNDaughters();
2078 for(Int_t i=0;i< nDaugthers;++i){
2079 AliAODMCParticle *partdaughter = (AliAODMCParticle*)(tca->At(i));
2080 if(!partdaughter) continue;
2081 if(partdaughter->Pt()<fTrackPtCut)continue;
2082 if(TMath::Abs(partdaughter->Eta())>fTrackEtaWindow)continue;
2083 if(partdaughter->Charge()==0)continue;
2085 if(!IsDMeson(partdaughter->PdgCode()) && !IsBMeson(partdaughter->PdgCode()) ){
2086 list->Add(partdaughter);
2089 else count+=AddDaughters(list,part,tca);
2096 void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() {
2099 AliInfo("Trying to connect to AliEn ...");
2100 TGrid::Connect("alien://");
2103 TFile *f = TFile::Open(fPathTrPtResolution.Data());
2107 TProfile *fProfPtPtSigma1PtGlobSt = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
2108 TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
2109 TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
2111 SetSmearResolution(kTRUE);
2112 SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD);
2117 void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() {
2120 AliInfo("Trying to connect to AliEn ...");
2121 TGrid::Connect("alien://");
2124 TFile *f = TFile::Open(fPathTrEfficiency.Data());
2127 TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
2128 TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
2129 TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
2131 SetDiceEfficiency(kTRUE);
2133 if(fChangeEfficiencyFraction>0.) {
2135 TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin");
2137 for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) {
2138 Double_t content = hEffPosGlobSt->GetBinContent(bin);
2139 hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction);
2142 SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
2146 SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
2150 void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
2153 // set mom res profiles
2156 if(fMomResH1) delete fMomResH1;
2157 if(fMomResH2) delete fMomResH2;
2158 if(fMomResH3) delete fMomResH3;
2160 fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
2161 fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
2162 fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
2166 void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
2168 // set tracking efficiency histos
2171 fhEffH1 = (TH1*)h1->Clone("fhEffH1");
2172 fhEffH2 = (TH1*)h2->Clone("fhEffH2");
2173 fhEffH3 = (TH1*)h3->Clone("fhEffH3");
2176 Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
2179 // Get smearing on generated momentum
2182 //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
2184 TProfile *fMomRes = 0x0;
2185 if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
2186 if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
2187 if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
2194 Double_t smear = 0.;
2197 if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
2198 if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
2199 if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
2203 Int_t bin = fMomRes->FindBin(pt);
2205 smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
2209 if(fMomRes) delete fMomRes;
2214 void AliAnalysisTaskJetCluster::FitMomentumResolution() {
2216 // Fit linear function on momentum resolution at high pT
2219 if(!fMomResH1Fit && fMomResH1) {
2220 fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
2221 fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
2222 fMomResH1Fit ->SetRange(5.,100.);
2225 if(!fMomResH2Fit && fMomResH2) {
2226 fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
2227 fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
2228 fMomResH2Fit ->SetRange(5.,100.);
2231 if(!fMomResH3Fit && fMomResH3) {
2232 fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
2233 fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
2234 fMomResH3Fit ->SetRange(5.,100.);
2240 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
2241 for(int i = 0; i < particles.GetEntries(); i++){
2242 AliVParticle *vp = (AliVParticle*)particles.At(i);
2243 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
2244 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
2245 jInp.set_user_index(i);
2246 inputParticles.push_back(jInp);
2255 bool AliAnalysisTaskJetCluster::IsBMeson(int pc){
2256 int bPdG[] = {511,521,10511,10521,513,523,10513,10523,20513,20523,20513,20523,515,525,531,
2257 10531,533,10533,20533,535,541,10541,543,10543,20543,545,551,10551,100551,
2258 110551,200551,210551,553,10553,20553,30553,100553,110553,120553,130553,200553,210553,220553,
2259 300553,9000533,9010553,555,10555,20555,100555,110555,120555,200555,557,100557};
2260 for(int i=0;i< (int)(sizeof(bPdG)/sizeof(int));++i) if(abs(pc) == bPdG[i]) return true;
2263 bool AliAnalysisTaskJetCluster::IsDMeson(int pc){
2264 int bPdG[] = {411,421,10411,10421,413,423,10413,10423,20431,20423,415,
2265 425,431,10431,433,10433,20433,435,441,10441,100441,443,10443,20443,
2266 100443,30443,9000443,9010443,9020443,445,100445};
2267 for(int i=0;i< (int)(sizeof(bPdG)/sizeof(int));++i) if(abs(pc) == bPdG[i]) return true;