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),
124 fJetOutputMinPt(0.150),
125 fMaxTrackPtInJet(100.),
131 fStoreRhoLeadingTrackCorr(kFALSE),
133 fBackgroundBranch(""),
144 fUseTrPtResolutionSmearing(kFALSE),
145 fUseDiceEfficiency(kFALSE),
146 fDiceEfficiencyMinPt(-1.),
147 fUseTrPtResolutionFromOADB(kFALSE),
148 fUseTrEfficiencyFromOADB(kFALSE),
149 fPathTrPtResolution(""),
150 fPathTrEfficiency(""),
151 fChangeEfficiencyFraction(0.),
152 fEfficiencyFixed(1.),
154 fAlgorithm(fastjet::kt_algorithm),
155 fStrategy(fastjet::Best),
156 fRecombScheme(fastjet::BIpt_scheme),
157 fAreaType(fastjet::active_area),
159 fActiveAreaRepeats(1),
163 fTCARandomConesOut(0x0),
164 fTCARandomConesOutRan(0x0),
165 fAODJetBackgroundOut(0x0),
171 fh1PtHardTrials(0x0),
174 fh1NConstLeadingRec(0x0),
176 fh1PtJetsLeadingRecIn(0x0),
177 fh1PtJetConstRec(0x0),
178 fh1PtJetConstLeadingRec(0x0),
179 fh1PtTracksRecIn(0x0),
180 fh1PtTracksLeadingRecIn(0x0),
182 fh1NConstRecRan(0x0),
183 fh1PtJetsLeadingRecInRan(0x0),
184 fh1NConstLeadingRecRan(0x0),
185 fh1PtJetsRecInRan(0x0),
186 fh1PtTracksGenIn(0x0),
188 fh1CentralityPhySel(0x0),
190 fh1CentralitySelect(0x0),
195 fh2NRecTracksPt(0x0),
197 fh2NConstLeadingPt(0x0),
199 fh2LeadingJetPhiEta(0x0),
201 fh2LeadingJetEtaPt(0x0),
203 fh2LeadingTrackEtaPt(0x0),
204 fh2JetsLeadingPhiEta(0x0),
205 fh2JetsLeadingPhiPt(0x0),
206 fh2TracksLeadingPhiEta(0x0),
207 fh2TracksLeadingPhiPt(0x0),
208 fh2TracksLeadingJetPhiPt(0x0),
209 fh2JetsLeadingPhiPtW(0x0),
210 fh2TracksLeadingPhiPtW(0x0),
211 fh2TracksLeadingJetPhiPtW(0x0),
212 fh2NRecJetsPtRan(0x0),
214 fh2NConstLeadingPtRan(0x0),
219 fh2TracksLeadingJetPhiPtRan(0x0),
220 fh2TracksLeadingJetPhiPtWRan(0x0),
221 fh3CentvsRhoLeadingTrackPt(0x0),
222 fh3CentvsSigmaLeadingTrackPt(0x0),
223 fh3MultvsRhoLeadingTrackPt(0x0),
224 fh3MultvsSigmaLeadingTrackPt(0x0),
225 fh3CentvsRhoLeadingTrackPtQ1(0x0),
226 fh3CentvsRhoLeadingTrackPtQ2(0x0),
227 fh3CentvsRhoLeadingTrackPtQ3(0x0),
228 fh3CentvsRhoLeadingTrackPtQ4(0x0),
229 fh3CentvsSigmaLeadingTrackPtQ1(0x0),
230 fh3CentvsSigmaLeadingTrackPtQ2(0x0),
231 fh3CentvsSigmaLeadingTrackPtQ3(0x0),
232 fh3CentvsSigmaLeadingTrackPtQ4(0x0),
233 fh3MultvsRhoLeadingTrackPtQ1(0x0),
234 fh3MultvsRhoLeadingTrackPtQ2(0x0),
235 fh3MultvsRhoLeadingTrackPtQ3(0x0),
236 fh3MultvsRhoLeadingTrackPtQ4(0x0),
237 fh3MultvsSigmaLeadingTrackPtQ1(0x0),
238 fh3MultvsSigmaLeadingTrackPtQ2(0x0),
239 fh3MultvsSigmaLeadingTrackPtQ3(0x0),
240 fh3MultvsSigmaLeadingTrackPtQ4(0x0),
241 fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0),
242 fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0),
243 fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0),
244 fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0),
245 fh2PtGenPtSmeared(0x0),
247 fp1PtResolution(0x0),
254 for(int i = 0;i<3;i++){
255 fh1BiARandomCones[i] = 0;
256 fh1BiARandomConesRan[i] = 0;
258 for(int i = 0;i<kMaxCent;i++){
259 fh2JetsLeadingPhiPtC[i] = 0;
260 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
261 fh2TracksLeadingJetPhiPtC[i] = 0;
262 fh2TracksLeadingJetPhiPtWC[i] = 0;
266 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
267 AliAnalysisTaskSE(name),
271 fUseAODTrackInput(kFALSE),
272 fUseAODMCInput(kFALSE),
273 fUseBackgroundCalc(kFALSE),
274 fEventSelection(kFALSE),
275 fRequireVZEROAC(kFALSE),
276 fRequireTZEROvtx(kFALSE),
278 fFilterMaskBestPt(0),
281 fTrackTypeRec(kTrackUndef),
282 fTrackTypeGen(kTrackUndef),
284 fNSkipLeadingCone(0),
288 fTrackEtaWindow(0.9),
292 fJetOutputMinPt(0.150),
293 fMaxTrackPtInJet(100.),
299 fStoreRhoLeadingTrackCorr(kFALSE),
301 fBackgroundBranch(""),
312 fUseTrPtResolutionSmearing(kFALSE),
313 fUseDiceEfficiency(kFALSE),
314 fDiceEfficiencyMinPt(-1.),
315 fUseTrPtResolutionFromOADB(kFALSE),
316 fUseTrEfficiencyFromOADB(kFALSE),
317 fPathTrPtResolution(""),
318 fPathTrEfficiency(""),
319 fChangeEfficiencyFraction(0.),
320 fEfficiencyFixed(1.),
322 fAlgorithm(fastjet::kt_algorithm),
323 fStrategy(fastjet::Best),
324 fRecombScheme(fastjet::BIpt_scheme),
325 fAreaType(fastjet::active_area),
327 fActiveAreaRepeats(1),
331 fTCARandomConesOut(0x0),
332 fTCARandomConesOutRan(0x0),
333 fAODJetBackgroundOut(0x0),
339 fh1PtHardTrials(0x0),
342 fh1NConstLeadingRec(0x0),
344 fh1PtJetsLeadingRecIn(0x0),
345 fh1PtJetConstRec(0x0),
346 fh1PtJetConstLeadingRec(0x0),
347 fh1PtTracksRecIn(0x0),
348 fh1PtTracksLeadingRecIn(0x0),
350 fh1NConstRecRan(0x0),
351 fh1PtJetsLeadingRecInRan(0x0),
352 fh1NConstLeadingRecRan(0x0),
353 fh1PtJetsRecInRan(0x0),
354 fh1PtTracksGenIn(0x0),
356 fh1CentralityPhySel(0x0),
358 fh1CentralitySelect(0x0),
363 fh2NRecTracksPt(0x0),
365 fh2NConstLeadingPt(0x0),
367 fh2LeadingJetPhiEta(0x0),
369 fh2LeadingJetEtaPt(0x0),
371 fh2LeadingTrackEtaPt(0x0),
372 fh2JetsLeadingPhiEta(0x0),
373 fh2JetsLeadingPhiPt(0x0),
374 fh2TracksLeadingPhiEta(0x0),
375 fh2TracksLeadingPhiPt(0x0),
376 fh2TracksLeadingJetPhiPt(0x0),
377 fh2JetsLeadingPhiPtW(0x0),
378 fh2TracksLeadingPhiPtW(0x0),
379 fh2TracksLeadingJetPhiPtW(0x0),
380 fh2NRecJetsPtRan(0x0),
382 fh2NConstLeadingPtRan(0x0),
387 fh2TracksLeadingJetPhiPtRan(0x0),
388 fh2TracksLeadingJetPhiPtWRan(0x0),
389 fh3CentvsRhoLeadingTrackPt(0x0),
390 fh3CentvsSigmaLeadingTrackPt(0x0),
391 fh3MultvsRhoLeadingTrackPt(0x0),
392 fh3MultvsSigmaLeadingTrackPt(0x0),
393 fh3CentvsRhoLeadingTrackPtQ1(0x0),
394 fh3CentvsRhoLeadingTrackPtQ2(0x0),
395 fh3CentvsRhoLeadingTrackPtQ3(0x0),
396 fh3CentvsRhoLeadingTrackPtQ4(0x0),
397 fh3CentvsSigmaLeadingTrackPtQ1(0x0),
398 fh3CentvsSigmaLeadingTrackPtQ2(0x0),
399 fh3CentvsSigmaLeadingTrackPtQ3(0x0),
400 fh3CentvsSigmaLeadingTrackPtQ4(0x0),
401 fh3MultvsRhoLeadingTrackPtQ1(0x0),
402 fh3MultvsRhoLeadingTrackPtQ2(0x0),
403 fh3MultvsRhoLeadingTrackPtQ3(0x0),
404 fh3MultvsRhoLeadingTrackPtQ4(0x0),
405 fh3MultvsSigmaLeadingTrackPtQ1(0x0),
406 fh3MultvsSigmaLeadingTrackPtQ2(0x0),
407 fh3MultvsSigmaLeadingTrackPtQ3(0x0),
408 fh3MultvsSigmaLeadingTrackPtQ4(0x0),
409 fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0),
410 fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0),
411 fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0),
412 fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0),
413 fh2PtGenPtSmeared(0x0),
415 fp1PtResolution(0x0),
422 for(int i = 0;i<3;i++){
423 fh1BiARandomCones[i] = 0;
424 fh1BiARandomConesRan[i] = 0;
426 for(int i = 0;i<kMaxCent;i++){
427 fh2JetsLeadingPhiPtC[i] = 0;
428 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
429 fh2TracksLeadingJetPhiPtC[i] = 0;
430 fh2TracksLeadingJetPhiPtWC[i] = 0;
432 DefineOutput(1, TList::Class());
437 Bool_t AliAnalysisTaskJetCluster::Notify()
440 // Implemented Notify() to read the cross sections
441 // and number of trials from pyxsec.root
446 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
450 // Create the output container
453 fRandom = new TRandom3(0);
459 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
463 if(fNonStdBranch.Length()!=0)
465 // only create the output branch if we have a name
466 // Create a new branch for jets...
467 // -> cleared in the UserExec....
468 // here we can also have the case that the brnaches are written to a separate file
471 fTCAJetsOut = new TClonesArray("AliAODJet", 0);
472 fTCAJetsOut->SetName(fNonStdBranch.Data());
473 AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
476 if(fJetTypes&kJetRan){
477 fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
478 fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
479 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
480 if( fEfficiencyFixed < 1.)
481 fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dEffFixed%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.)));
483 fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dFr%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
485 AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
488 if(fUseBackgroundCalc){
489 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
490 fAODJetBackgroundOut = new AliAODJetEventBackground();
491 fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
492 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
493 if( fEfficiencyFixed < 1.)
494 fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dEffFixed%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.)));
496 fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dFr%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
498 AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
501 // create the branch for the random cones with the same R
502 TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
503 if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
504 if( fEfficiencyFixed < 1.)
505 cName = Form("%sDetector%d%dEffFixed%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.),fNSkipLeadingCone);
507 cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone);
511 if(!AODEvent()->FindListObject(cName.Data())){
512 fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
513 fTCARandomConesOut->SetName(cName.Data());
514 AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
517 // create the branch with the random for the random cones on the random event
518 if(fJetTypes&kRCRan){
519 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
520 if(!AODEvent()->FindListObject(cName.Data())){
521 fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
522 fTCARandomConesOutRan->SetName(cName.Data());
523 AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
528 if(fNonStdFile.Length()!=0){
530 // case that we have an AOD extension we need to fetch the jets from the extended output
531 // we identify the extension aod event by looking for the branchname
532 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
533 // case that we have an AOD extension we need can fetch the background maybe from the extended output
534 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
539 if(!fHistList)fHistList = new TList();
540 fHistList->SetOwner();
541 PostData(1, fHistList); // post data in any case once
543 Bool_t oldStatus = TH1::AddDirectoryStatus();
544 TH1::AddDirectory(kFALSE);
549 const Int_t nBinPt = 100;
550 Double_t binLimitsPt[nBinPt+1];
551 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
553 binLimitsPt[iPt] = 0.0;
556 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
560 const Int_t nBinPhi = 90;
561 Double_t binLimitsPhi[nBinPhi+1];
562 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
564 binLimitsPhi[iPhi] = -1.*TMath::Pi();
567 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
573 const Int_t nBinEta = 40;
574 Double_t binLimitsEta[nBinEta+1];
575 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
577 binLimitsEta[iEta] = -2.0;
580 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
584 const int nChMax = 5000;
586 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
587 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
589 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
590 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
593 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
594 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
596 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
597 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
598 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
599 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
602 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
603 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
604 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
606 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
607 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
608 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
609 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
610 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
611 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
612 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
613 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
614 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
615 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
617 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
618 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
619 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
621 fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
622 fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
623 fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
625 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
626 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
627 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
631 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
632 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
633 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
634 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
636 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
637 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);
638 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);
639 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);
643 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
644 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
645 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
646 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
648 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
649 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
650 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
651 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
653 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
654 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
655 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
656 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
660 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
661 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
662 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
663 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
664 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
665 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
666 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
667 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
668 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
669 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
670 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
671 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
673 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
674 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
675 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
676 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
678 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
679 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
680 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
681 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
683 if(fStoreRhoLeadingTrackCorr) {
684 fh3CentvsRhoLeadingTrackPt = new TH3F("fh3CentvsRhoLeadingTrackPt","centrality vs background density full event; centrality; #rho", 50,0.,100.,500,0.,250.,100,0.,100.);
685 fh3CentvsSigmaLeadingTrackPt = new TH3F("fh3CentvsSigmaLeadingTrackPt","centrality vs sigma full event; centrality; #sigma(#rho)", 50,0.,100.,50,0.,50.,100,0.,100.);
686 fh3MultvsRhoLeadingTrackPt = new TH3F("fh3MultvsRhoLeadingTrackPt","multiplicity vs background density full event; multiplicity; #rho", 100,0.,5000.,500,0.,250.,100,0.,100.);
687 fh3MultvsSigmaLeadingTrackPt = new TH3F("fh3MultvsSigmaLeadingTrackPt","multiplicity vs sigma full event; multiplicity; #sigma(#rho)", 100,0.,5000.,50,0.,50.,100,0.,100.);
690 fh3CentvsRhoLeadingTrackPtQ1 = new TH3F("fh3CentvsRhoLeadingTrackPtQ1","centrality vs background density Q1; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
691 fh3CentvsRhoLeadingTrackPtQ2 = new TH3F("fh3CentvsRhoLeadingTrackPtQ2","centrality vs background density Q2; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
692 fh3CentvsRhoLeadingTrackPtQ3 = new TH3F("fh3CentvsRhoLeadingTrackPtQ3","centrality vs background density Q3; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
693 fh3CentvsRhoLeadingTrackPtQ4 = new TH3F("fh3CentvsRhoLeadingTrackPtQ4","centrality vs background density Q4; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
695 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.);
696 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.);
697 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.);
698 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.);
700 fh3MultvsRhoLeadingTrackPtQ1 = new TH3F("fh3MultvsRhoLeadingTrackPtQ1","multiplicity vs background density Q1; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
701 fh3MultvsRhoLeadingTrackPtQ2 = new TH3F("fh3MultvsRhoLeadingTrackPtQ2","multiplicity vs background density Q2; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
702 fh3MultvsRhoLeadingTrackPtQ3 = new TH3F("fh3MultvsRhoLeadingTrackPtQ3","multiplicity vs background density Q3; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
703 fh3MultvsRhoLeadingTrackPtQ4 = new TH3F("fh3MultvsRhoLeadingTrackPtQ4","multiplicity vs background density Q4; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
705 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.);
706 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.);
707 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.);
708 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.);
711 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.);
712 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.);
713 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.);
714 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.);
716 fHistList->Add(fh3CentvsRhoLeadingTrackPt);
717 fHistList->Add(fh3CentvsSigmaLeadingTrackPt);
718 fHistList->Add(fh3MultvsRhoLeadingTrackPt);
719 fHistList->Add(fh3MultvsSigmaLeadingTrackPt);
721 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ1);
722 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ2);
723 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ3);
724 fHistList->Add(fh3CentvsRhoLeadingTrackPtQ4);
726 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ1);
727 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ2);
728 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ3);
729 fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ4);
731 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ1);
732 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ2);
733 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ3);
734 fHistList->Add(fh3MultvsRhoLeadingTrackPtQ4);
736 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ1);
737 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ2);
738 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ3);
739 fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ4);
741 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ1);
742 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ2);
743 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ3);
744 fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ4);
748 //Detector level effects histos
749 fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
751 fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
752 fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
754 fHistList->Add(fh2PtGenPtSmeared);
755 fHistList->Add(fp1Efficiency);
756 fHistList->Add(fp1PtResolution);
758 if(fNRandomCones>0&&fUseBackgroundCalc){
759 for(int i = 0;i<3;i++){
760 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
761 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
765 for(int i = 0;i < kMaxCent;i++){
766 fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
767 fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
768 fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
769 fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
772 const Int_t saveLevel = 3; // large save level more histos
774 fHistList->Add(fh1Xsec);
775 fHistList->Add(fh1Trials);
777 fHistList->Add(fh1NJetsRec);
778 fHistList->Add(fh1NConstRec);
779 fHistList->Add(fh1NConstLeadingRec);
780 fHistList->Add(fh1PtJetsRecIn);
781 fHistList->Add(fh1PtJetsLeadingRecIn);
782 fHistList->Add(fh1PtTracksRecIn);
783 fHistList->Add(fh1PtTracksLeadingRecIn);
784 fHistList->Add(fh1PtJetConstRec);
785 fHistList->Add(fh1PtJetConstLeadingRec);
786 fHistList->Add(fh1NJetsRecRan);
787 fHistList->Add(fh1NConstRecRan);
788 fHistList->Add(fh1PtJetsLeadingRecInRan);
789 fHistList->Add(fh1NConstLeadingRecRan);
790 fHistList->Add(fh1PtJetsRecInRan);
791 fHistList->Add(fh1Nch);
792 fHistList->Add(fh1Centrality);
793 fHistList->Add(fh1CentralitySelect);
794 fHistList->Add(fh1CentralityPhySel);
795 fHistList->Add(fh1Z);
796 fHistList->Add(fh1ZSelect);
797 fHistList->Add(fh1ZPhySel);
798 if(fNRandomCones>0&&fUseBackgroundCalc){
799 for(int i = 0;i<3;i++){
800 fHistList->Add(fh1BiARandomCones[i]);
801 fHistList->Add(fh1BiARandomConesRan[i]);
804 for(int i = 0;i < kMaxCent;i++){
805 fHistList->Add(fh2JetsLeadingPhiPtC[i]);
806 fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
807 fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
808 fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
811 fHistList->Add(fh2NRecJetsPt);
812 fHistList->Add(fh2NRecTracksPt);
813 fHistList->Add(fh2NConstPt);
814 fHistList->Add(fh2NConstLeadingPt);
815 fHistList->Add(fh2PtNch);
816 fHistList->Add(fh2PtNchRan);
817 fHistList->Add(fh2PtNchN);
818 fHistList->Add(fh2PtNchNRan);
819 fHistList->Add(fh2JetPhiEta);
820 fHistList->Add(fh2LeadingJetPhiEta);
821 fHistList->Add(fh2JetEtaPt);
822 fHistList->Add(fh2LeadingJetEtaPt);
823 fHistList->Add(fh2TrackEtaPt);
824 fHistList->Add(fh2LeadingTrackEtaPt);
825 fHistList->Add(fh2JetsLeadingPhiEta );
826 fHistList->Add(fh2JetsLeadingPhiPt);
827 fHistList->Add(fh2TracksLeadingPhiEta);
828 fHistList->Add(fh2TracksLeadingPhiPt);
829 fHistList->Add(fh2TracksLeadingJetPhiPt);
830 fHistList->Add(fh2JetsLeadingPhiPtW);
831 fHistList->Add(fh2TracksLeadingPhiPtW);
832 fHistList->Add(fh2TracksLeadingJetPhiPtW);
833 fHistList->Add(fh2NRecJetsPtRan);
834 fHistList->Add(fh2NConstPtRan);
835 fHistList->Add(fh2NConstLeadingPtRan);
836 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
837 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
840 // =========== Switch on Sumw2 for all histos ===========
841 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
842 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
847 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
850 TH1::AddDirectory(oldStatus);
853 void AliAnalysisTaskJetCluster::LocalInit()
859 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
861 if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB();
862 if(fUseTrEfficiencyFromOADB) LoadTrEfficiencyRootFileFromOADB();
865 FitMomentumResolution();
869 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
872 // handle and reset the output jet branch
874 if(fTCAJetsOut)fTCAJetsOut->Delete();
875 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
876 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
877 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
878 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
880 AliAODJetEventBackground* externalBackground = 0;
881 if(!externalBackground&&fBackgroundBranch.Length()){
882 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
883 if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
884 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
887 // Execute analysis for current event
889 AliESDEvent *fESD = 0;
890 if(fUseAODTrackInput){
891 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
893 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
899 // assume that the AOD is in the general output...
902 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
906 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
910 //Check if information is provided detector level effects
911 if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE;
912 if(!fhEffH1 || !fhEffH2 || !fhEffH3 ) fUseDiceEfficiency = kFALSE;
913 if( fEfficiencyFixed < 1. ) fUseDiceEfficiency = kTRUE;
915 Bool_t selectEvent = false;
916 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
922 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
923 TString vtxTitle(vtxAOD->GetTitle());
924 zVtx = vtxAOD->GetZ();
926 cent = fAOD->GetHeader()->GetCentrality();
927 if(cent<10)cenClass = 0;
928 else if(cent<30)cenClass = 1;
929 else if(cent<50)cenClass = 2;
930 else if(cent<80)cenClass = 3;
931 if(physicsSelection){
932 fh1CentralityPhySel->Fill(cent);
933 fh1ZPhySel->Fill(zVtx);
937 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
938 Float_t yvtx = vtxAOD->GetY();
939 Float_t xvtx = vtxAOD->GetX();
940 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
941 if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
942 if(physicsSelection){
948 if(cent<fCentCutLo||cent>fCentCutUp){
960 const AliAODVZERO *vzero = fAOD->GetVZEROData();
962 if((vzero->GetTriggerChargeA()>0)&&(vzero->GetTriggerChargeC()>0)){
967 const AliAODTZERO *tzero = fAOD->GetTZEROData();
969 if(TMath::Abs(tzero->GetT0VertexRaw())<100){
974 if(fRequireVZEROAC&&fRequireTZEROvtx)selectEvent = selectEvent&&V0&&T0;
975 else if(fRequireTZEROvtx)selectEvent = selectEvent&&T0;
976 else if(fRequireVZEROAC)selectEvent = selectEvent&&V0;
980 PostData(1, fHistList);
983 fh1Centrality->Fill(cent);
985 fh1Trials->Fill("#sum{ntrials}",1);
988 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
990 // ==== General variables needed
994 // we simply fetch the tracks/mc particles as a list of AliVParticles
999 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
1000 Float_t nCh = recParticles.GetEntries();
1002 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
1003 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
1004 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
1006 // find the jets....
1008 vector<fastjet::PseudoJet> inputParticlesRec;
1009 vector<fastjet::PseudoJet> inputParticlesRecRan;
1011 // Generate the random cones
1013 AliAODJet vTmpRan(1,0,0,1);
1014 for(int i = 0; i < recParticles.GetEntries(); i++){
1015 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1017 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1018 // we take total momentum here
1020 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
1021 //Add particles to fastjet in case we are not running toy model
1022 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
1023 jInp.set_user_index(i);
1024 inputParticlesRec.push_back(jInp);
1026 else if(fUseDiceEfficiency) {
1028 // Dice to decide if particle is kept or not - toy model for efficiency
1030 Double_t sumEff = 0.;
1032 Double_t eff[3] = {0.};
1035 Double_t rnd = fRandom->Uniform(1.);
1036 if( fEfficiencyFixed<1. ) {
1037 sumEff = fEfficiencyFixed;
1041 Double_t pTtmp = pT;
1042 if(pT>10.) pTtmp = 10.;
1043 Double_t eff1 = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
1044 Double_t eff2 = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
1045 Double_t eff3 = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));
1047 //Sort efficiencies from large to small
1048 if(eff1>eff2 && eff1>eff3) {
1063 else if(eff2>eff1 && eff2>eff3) {
1078 else if(eff3>eff1 && eff3>eff2) {
1094 sumEff = eff[0]+eff[1]+eff[2];
1096 fp1Efficiency->Fill(vp->Pt(),sumEff);
1097 if(rnd>sumEff && pT > fDiceEfficiencyMinPt) continue;
1099 if(fUseTrPtResolutionSmearing) {
1100 //Smear momentum of generated particle
1101 Double_t smear = 1.;
1102 //Select hybrid track category
1104 smear = GetMomentumSmearing(cat[2],pT);
1105 else if(rnd<=(eff[2]+eff[1]))
1106 smear = GetMomentumSmearing(cat[1],pT);
1108 smear = GetMomentumSmearing(cat[0],pT);
1110 fp1PtResolution->Fill(vp->Pt(),smear);
1112 Double_t sigma = vp->Pt()*smear;
1113 Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
1115 Double_t phi = vp->Phi();
1116 Double_t theta = 2.*TMath::ATan(TMath::Exp(-1.*(vp->Eta())));
1117 Double_t pX = pTrec * TMath::Cos(phi);
1118 Double_t pY = pTrec * TMath::Sin(phi);
1119 Double_t pZ = pTrec/TMath::Tan(theta);
1120 Double_t p=TMath::Sqrt(pTrec*pTrec+pZ*pZ);
1122 fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);
1124 fastjet::PseudoJet jInp(pX,pY,pZ,p);
1125 jInp.set_user_index(i);
1126 inputParticlesRec.push_back(jInp);
1130 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
1131 jInp.set_user_index(i);
1132 inputParticlesRec.push_back(jInp);
1138 // the randomized input changes eta and phi, but keeps the p_T
1139 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1140 Double_t pT = vp->Pt();
1141 Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
1142 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
1144 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
1145 Double_t pZ = pT/TMath::Tan(theta);
1147 Double_t pX = pT * TMath::Cos(phi);
1148 Double_t pY = pT * TMath::Sin(phi);
1149 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
1150 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
1152 jInpRan.set_user_index(i);
1153 inputParticlesRecRan.push_back(jInpRan);
1154 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
1157 // fill the tref array, only needed when we write out jets
1160 fRef->Delete(); // make sure to delete before placement new...
1161 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
1162 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
1165 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ...
1169 if(inputParticlesRec.size()==0){
1170 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
1171 PostData(1, fHistList);
1176 // employ setters for these...
1179 // now create the object that holds info about ghosts
1181 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
1182 // reduce CPU time...
1184 fActiveAreaRepeats = 0;
1188 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
1189 fastjet::AreaType areaType = fastjet::active_area;
1190 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
1191 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
1192 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
1194 //range where to compute background
1195 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
1197 phiMax = 2*TMath::Pi();
1198 rapMax = fGhostEtamax - fRparam;
1199 rapMin = - fGhostEtamax + fRparam;
1200 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
1203 const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
1204 const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
1207 fh1NJetsRec->Fill(sortedJets.size());
1209 // loop over all jets an fill information, first on is the leading jet
1211 Int_t nRecOver = inclusiveJets.size();
1212 Int_t nRec = inclusiveJets.size();
1213 if(inclusiveJets.size()>0){
1214 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
1215 Double_t area = clustSeq.area(sortedJets[0]);
1216 leadingJet.SetEffArea(area,0);
1217 Float_t pt = leadingJet.Pt();
1218 Int_t nAodOutJets = 0;
1219 Int_t nAodOutTracks = 0;
1220 AliAODJet *aodOutJet = 0;
1223 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
1224 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
1225 while(pt<ptCut&&iCount<nRec){
1229 pt = sortedJets[iCount].perp();
1232 if(nRecOver<=0)break;
1233 fh2NRecJetsPt->Fill(ptCut,nRecOver);
1235 Float_t phi = leadingJet.Phi();
1236 if(phi<0)phi+=TMath::Pi()*2.;
1237 Float_t eta = leadingJet.Eta();
1239 if(externalBackground){
1240 // carefull has to be filled in a task before
1241 // todo, ReArrange to the botom
1242 pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
1244 pt = leadingJet.Pt() - pTback;
1245 // correlation of leading jet with tracks
1246 TIterator *recIter = recParticles.MakeIterator();
1248 AliVParticle *tmpRecTrack = 0;
1249 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
1250 Float_t tmpPt = tmpRecTrack->Pt();
1252 Float_t tmpPhi = tmpRecTrack->Phi();
1253 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1254 Float_t dPhi = phi - tmpPhi;
1255 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1256 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1257 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
1258 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
1260 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
1261 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1267 TLorentzVector vecareab;
1268 for(int j = 0; j < nRec;j++){
1269 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
1272 Float_t tmpPt = tmpRec.Pt();
1274 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
1275 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
1276 aodOutJet->GetRefTracks()->Clear();
1277 Double_t area1 = clustSeq.area(sortedJets[j]);
1278 aodOutJet->SetEffArea(area1,0);
1279 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
1280 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
1281 aodOutJet->SetVectorAreaCharged(&vecareab);
1285 Float_t tmpPtBack = 0;
1286 if(externalBackground){
1287 // carefull has to be filled in a task before
1288 // todo, ReArrange to the botom
1289 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
1291 tmpPt = tmpPt - tmpPtBack;
1292 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
1294 fh1PtJetsRecIn->Fill(tmpPt);
1295 // Fill Spectra with constituentsemacs
1296 const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
1298 fh1NConstRec->Fill(constituents.size());
1299 fh2PtNch->Fill(nCh,tmpPt);
1300 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
1301 fh2NConstPt->Fill(tmpPt,constituents.size());
1302 // loop over constiutents and fill spectrum
1304 AliVParticle *partLead = 0;
1305 Float_t ptLead = -1;
1307 for(unsigned int ic = 0; ic < constituents.size();ic++){
1308 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1310 fh1PtJetConstRec->Fill(part->Pt());
1312 if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1313 if(part->Pt()>fMaxTrackPtInJet){
1314 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1317 if(part->Pt()>ptLead){
1318 ptLead = part->Pt();
1321 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1324 AliAODTrack *aodT = 0;
1327 //set pT of leading constituent of jet
1328 aodOutJet->SetPtLeading(partLead->Pt());
1329 aodT = dynamic_cast<AliAODTrack*>(partLead);
1331 if(aodT->TestFilterBit(fFilterMaskBestPt)){
1332 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
1339 Float_t tmpPhi = tmpRec.Phi();
1340 Float_t tmpEta = tmpRec.Eta();
1341 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1343 fh1PtJetsLeadingRecIn->Fill(tmpPt);
1344 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1345 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1346 fh1NConstLeadingRec->Fill(constituents.size());
1347 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1350 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1351 fh2JetEtaPt->Fill(tmpEta,tmpPt);
1352 Float_t dPhi = phi - tmpPhi;
1353 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1354 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1355 Float_t dEta = eta - tmpRec.Eta();
1356 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1357 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1359 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1360 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1362 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1363 }// loop over reconstructed jets
1368 // Add the random cones...
1369 if(fNRandomCones>0&&fTCARandomConesOut){
1370 // create a random jet within the acceptance
1371 Double_t etaMax = fTrackEtaWindow - fRparam;
1374 Double_t pTC = 1; // small number
1375 for(int ir = 0;ir < fNRandomCones;ir++){
1376 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1377 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1379 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1380 Double_t pZC = pTC/TMath::Tan(thetaC);
1381 Double_t pXC = pTC * TMath::Cos(phiC);
1382 Double_t pYC = pTC * TMath::Sin(phiC);
1383 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1384 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
1386 for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1387 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1388 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1393 // test for overlap with previous cones to avoid double counting
1394 for(int iic = 0;iic<ir;iic++){
1395 AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1397 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1404 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1405 tmpRecC.SetPtLeading(-1.);
1406 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1407 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1408 }// loop over random cones creation
1411 // loop over the reconstructed particles and add up the pT in the random cones
1412 // maybe better to loop over randomized particles not in the real jets...
1413 // but this by definition brings dow average energy in the whole event
1414 AliAODJet vTmpRanR(1,0,0,1);
1415 for(int i = 0; i < recParticles.GetEntries(); i++){
1416 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1417 if(fTCARandomConesOut){
1418 for(int ir = 0;ir < fNRandomCones;ir++){
1419 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1420 if(jC&&jC->DeltaR(vp)<fRparam){
1421 if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1422 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1423 if(vp->Pt() > jC->GetPtLeading()) jC->SetPtLeading(vp->Pt());
1426 }// add up energy in cone
1428 // the randomized input changes eta and phi, but keeps the p_T
1429 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1430 Double_t pTR = vp->Pt();
1431 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1432 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1434 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1435 Double_t pZR = pTR/TMath::Tan(thetaR);
1437 Double_t pXR = pTR * TMath::Cos(phiR);
1438 Double_t pYR = pTR * TMath::Sin(phiR);
1439 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1440 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1441 if(fTCARandomConesOutRan){
1442 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1443 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1444 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1445 if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1446 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1447 if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt());
1452 }// loop over recparticles
1454 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1455 if(fTCARandomConesOut){
1456 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1457 // rescale the momentum vectors for the random cones
1459 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1461 Double_t etaC = rC->Eta();
1462 Double_t phiC = rC->Phi();
1463 // massless jet, unit vector
1464 pTC = rC->ChargedBgEnergy();
1465 if(pTC<=0)pTC = 0.001; // for almost empty events
1466 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1467 Double_t pZC = pTC/TMath::Tan(thetaC);
1468 Double_t pXC = pTC * TMath::Cos(phiC);
1469 Double_t pYC = pTC * TMath::Sin(phiC);
1470 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1471 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1472 rC->SetBgEnergy(0,0);
1473 rC->SetEffArea(jetArea,0);
1477 if(fTCARandomConesOutRan){
1478 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1479 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1482 Double_t etaC = rC->Eta();
1483 Double_t phiC = rC->Phi();
1484 // massless jet, unit vector
1485 pTC = rC->ChargedBgEnergy();
1486 if(pTC<=0)pTC = 0.001;// for almost empty events
1487 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1488 Double_t pZC = pTC/TMath::Tan(thetaC);
1489 Double_t pXC = pTC * TMath::Cos(phiC);
1490 Double_t pYC = pTC * TMath::Sin(phiC);
1491 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1492 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1493 rC->SetBgEnergy(0,0);
1494 rC->SetEffArea(jetArea,0);
1498 }// if(fNRandomCones
1500 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1502 if(fAODJetBackgroundOut){
1503 vector<fastjet::PseudoJet> jets2=sortedJets;
1504 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1508 Double_t meanarea1=0.;
1511 Double_t meanarea2=0.;
1515 Double_t meanarea4=0.;
1517 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1518 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1520 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1521 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1522 clustSeq.get_median_rho_and_sigma(sortedJets, range, true, bkg4, sigma4, meanarea4, true);
1523 fAODJetBackgroundOut->SetBackground(3,bkg4,sigma4,meanarea4);
1525 //////////////////////
1527 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1528 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1529 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1530 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1535 if(fStoreRhoLeadingTrackCorr) {
1536 vector<fastjet::PseudoJet> jets3=sortedJets;
1537 if(jets3.size()>2) jets3.erase(jets3.begin(),jets3.begin()+2);
1541 Double_t meanarea2=0.;
1543 clustSeq.get_median_rho_and_sigma(jets3, range, false, bkg2, sigma2, meanarea2, true);
1544 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1546 //Get leading track in event
1547 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1549 fh3CentvsRhoLeadingTrackPt->Fill(cent,bkg2,leading->Pt());
1550 fh3CentvsSigmaLeadingTrackPt->Fill(cent,sigma2,leading->Pt());
1551 fh3MultvsRhoLeadingTrackPt->Fill(nCh,bkg2,leading->Pt());
1552 fh3MultvsSigmaLeadingTrackPt->Fill(nCh,sigma2,leading->Pt());
1555 //Calculate rho with 4-vector area method for different quadrants with respect to the leading track in the event
1556 //exclude 2 leading jets from event
1557 //Quadrant 1: |phi_leadingTrack - phi_bkgCluster| < pi/2./2. - R (Near side to leading track)
1558 //Quadrant 2: pi/2 - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < pi/2 + (pi/2./2. - R) (Orthogonal to leading track)
1559 //Quadrant 3: pi - (pi/2/2 - R) < |phi_leadingTrack - phi_bkgCluster| < pi + (pi/2./2. - R) (Away side to leading track)
1560 //Quadrant 4: 3/2*pi - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < 3/2*pi + (pi/2./2. - R) (Orthogonal to leading track)
1562 //Assuming R=0.2 for background clusters
1564 Double_t bkg2Q[4] = {0.};
1565 Double_t sigma2Q[4] = {0.};
1566 Double_t meanarea2Q[4] = {0.};
1568 //Get phi of leading track in event
1569 Float_t phiLeadingTrack = leading->Phi();
1570 Float_t phiStep = (TMath::Pi()/2./2. - 0.2);
1572 //Quadrant1 - near side
1573 phiMin = phiLeadingTrack - phiStep;
1574 phiMax = phiLeadingTrack + phiStep;
1575 fastjet::RangeDefinition rangeQ1(rapMin,rapMax, phiMin, phiMax);
1576 clustSeq.get_median_rho_and_sigma(jets3, rangeQ1, false, bkg2Q[0], sigma2Q[0], meanarea2Q[0], true);
1578 //Quadrant2 - orthogonal
1579 Double_t phiQ2 = phiLeadingTrack + TMath::Pi()/2.;
1580 if(phiQ2>TMath::TwoPi()) phiQ2 = phiQ2 - TMath::TwoPi();
1581 phiMin = phiQ2 - phiStep;
1582 phiMax = phiQ2 + phiStep;
1583 fastjet::RangeDefinition rangeQ2(rapMin,rapMax, phiMin, phiMax);
1584 clustSeq.get_median_rho_and_sigma(jets3, rangeQ2, false, bkg2Q[1], sigma2Q[1], meanarea2Q[1], true);
1586 //Quadrant3 - away side
1587 Double_t phiQ3 = phiLeadingTrack + TMath::Pi();
1588 if(phiQ3>TMath::TwoPi()) phiQ3 = phiQ3 - TMath::TwoPi();
1589 phiMin = phiQ3 - phiStep;
1590 phiMax = phiQ3 + phiStep;
1591 fastjet::RangeDefinition rangeQ3(rapMin,rapMax, phiMin, phiMax);
1592 clustSeq.get_median_rho_and_sigma(jets3, rangeQ3, false, bkg2Q[2], sigma2Q[2], meanarea2Q[2], true);
1594 //Quadrant4 - othogonal
1595 Double_t phiQ4 = phiLeadingTrack + 3./2.*TMath::Pi();
1596 if(phiQ4>TMath::TwoPi()) phiQ4 = phiQ4 - TMath::TwoPi();
1597 phiMin = phiQ4 - phiStep;
1598 phiMax = phiQ4 + phiStep;
1599 fastjet::RangeDefinition rangeQ4(rapMin,rapMax, phiMin, phiMax);
1600 clustSeq.get_median_rho_and_sigma(jets3, rangeQ4, false, bkg2Q[3], sigma2Q[3], meanarea2Q[3], true);
1602 fh3CentvsRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0],leading->Pt());
1603 fh3CentvsRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1],leading->Pt());
1604 fh3CentvsRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2],leading->Pt());
1605 fh3CentvsRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3],leading->Pt());
1607 fh3CentvsSigmaLeadingTrackPtQ1->Fill(cent,sigma2Q[0],leading->Pt());
1608 fh3CentvsSigmaLeadingTrackPtQ2->Fill(cent,sigma2Q[1],leading->Pt());
1609 fh3CentvsSigmaLeadingTrackPtQ3->Fill(cent,sigma2Q[2],leading->Pt());
1610 fh3CentvsSigmaLeadingTrackPtQ4->Fill(cent,sigma2Q[3],leading->Pt());
1612 fh3MultvsRhoLeadingTrackPtQ1->Fill(nCh,bkg2Q[0],leading->Pt());
1613 fh3MultvsRhoLeadingTrackPtQ2->Fill(nCh,bkg2Q[1],leading->Pt());
1614 fh3MultvsRhoLeadingTrackPtQ3->Fill(nCh,bkg2Q[2],leading->Pt());
1615 fh3MultvsRhoLeadingTrackPtQ4->Fill(nCh,bkg2Q[3],leading->Pt());
1617 fh3MultvsSigmaLeadingTrackPtQ1->Fill(nCh,sigma2Q[0],leading->Pt());
1618 fh3MultvsSigmaLeadingTrackPtQ2->Fill(nCh,sigma2Q[1],leading->Pt());
1619 fh3MultvsSigmaLeadingTrackPtQ3->Fill(nCh,sigma2Q[2],leading->Pt());
1620 fh3MultvsSigmaLeadingTrackPtQ4->Fill(nCh,sigma2Q[3],leading->Pt());
1622 fh3CentvsDeltaRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0]-bkg2,leading->Pt());
1623 fh3CentvsDeltaRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1]-bkg2,leading->Pt());
1624 fh3CentvsDeltaRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2]-bkg2,leading->Pt());
1625 fh3CentvsDeltaRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3]-bkg2,leading->Pt());
1633 // fill track information
1634 Int_t nTrackOver = recParticles.GetSize();
1635 // do the same for tracks and jets
1638 TIterator *recIter = recParticles.MakeIterator();
1639 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1640 Float_t pt = tmpRec->Pt();
1642 // Printf("Leading track p_t %3.3E",pt);
1643 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1644 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1645 while(pt<ptCut&&tmpRec){
1647 tmpRec = (AliVParticle*)(recIter->Next());
1652 if(nTrackOver<=0)break;
1653 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1657 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1658 Float_t phi = leading->Phi();
1659 if(phi<0)phi+=TMath::Pi()*2.;
1660 Float_t eta = leading->Eta();
1662 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1663 Float_t tmpPt = tmpRec->Pt();
1664 Float_t tmpEta = tmpRec->Eta();
1665 fh1PtTracksRecIn->Fill(tmpPt);
1666 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1667 if(tmpRec==leading){
1668 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1669 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1673 Float_t tmpPhi = tmpRec->Phi();
1675 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1676 Float_t dPhi = phi - tmpPhi;
1677 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1678 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1679 Float_t dEta = eta - tmpRec->Eta();
1680 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1681 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1682 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1687 // find the random jets
1689 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1691 // fill the jet information from random track
1692 const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1693 const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1695 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1697 // loop over all jets an fill information, first on is the leading jet
1699 Int_t nRecOverRan = inclusiveJetsRan.size();
1700 Int_t nRecRan = inclusiveJetsRan.size();
1702 if(inclusiveJetsRan.size()>0){
1703 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1704 Float_t pt = leadingJet.Pt();
1707 TLorentzVector vecarearanb;
1709 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1710 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1711 while(pt<ptCut&&iCount<nRecRan){
1715 pt = sortedJetsRan[iCount].perp();
1718 if(nRecOverRan<=0)break;
1719 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1721 Float_t phi = leadingJet.Phi();
1722 if(phi<0)phi+=TMath::Pi()*2.;
1723 pt = leadingJet.Pt();
1725 // correlation of leading jet with random tracks
1727 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1729 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1731 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1732 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1733 Float_t dPhi = phi - tmpPhi;
1734 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1735 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1736 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1737 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1740 Int_t nAodOutJetsRan = 0;
1741 AliAODJet *aodOutJetRan = 0;
1742 for(int j = 0; j < nRecRan;j++){
1743 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1744 Float_t tmpPt = tmpRec.Pt();
1745 fh1PtJetsRecInRan->Fill(tmpPt);
1746 // Fill Spectra with constituents
1747 const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1748 fh1NConstRecRan->Fill(constituents.size());
1749 fh2NConstPtRan->Fill(tmpPt,constituents.size());
1750 fh2PtNchRan->Fill(nCh,tmpPt);
1751 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1754 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1755 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1756 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1757 aodOutJetRan->GetRefTracks()->Clear();
1758 aodOutJetRan->SetEffArea(arearan,0);
1759 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1760 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1761 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1766 Float_t tmpPhi = tmpRec.Phi();
1767 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1770 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1771 fh1NConstLeadingRecRan->Fill(constituents.size());
1772 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1778 if(fAODJetBackgroundOut){
1781 Double_t meanarea3=0.;
1782 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1783 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1784 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1786 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1787 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1796 // do the event selection if activated
1797 if(fJetTriggerPtCut>0){
1798 bool select = false;
1799 Float_t minPt = fJetTriggerPtCut;
1801 // hard coded for now ...
1802 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1803 if(cent<10)minPt = 50;
1804 else if(cent<30)minPt = 42;
1805 else if(cent<50)minPt = 28;
1806 else if(cent<80)minPt = 18;
1809 if(externalBackground)rho = externalBackground->GetBackground(2);
1811 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1812 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1813 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1822 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1823 fh1CentralitySelect->Fill(cent);
1824 fh1ZSelect->Fill(zVtx);
1825 aodH->SetFillAOD(kTRUE);
1829 if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1830 if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1831 if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1832 if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1834 PostData(1, fHistList);
1837 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1840 // Terminate analysis
1842 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1844 if(fMomResH1Fit) delete fMomResH1Fit;
1845 if(fMomResH2Fit) delete fMomResH2Fit;
1846 if(fMomResH3Fit) delete fMomResH3Fit;
1851 Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1854 // get list of tracks/particles for different types
1857 if(fDebug>2) Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1860 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly || type==kTrackAODMCextra || type==kTrackAODMCextraonly){
1862 if(type!=kTrackAODextraonly && type!=kTrackAODMCextraonly) {
1863 AliAODEvent *aod = 0;
1864 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1865 else aod = AODEvent();
1867 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1871 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1872 AliAODTrack *tr = aod->GetTrack(it);
1873 Bool_t bGood = false;
1874 if(fFilterType == 0)bGood = true;
1875 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1876 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1877 if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1878 if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1881 if(fRequireITSRefit==0){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
1882 if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1883 if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1886 if(tr->Pt()<fTrackPtCut){
1887 if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1890 if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1895 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1896 AliAODEvent *aod = 0;
1897 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1898 else aod = AODEvent();
1903 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1904 if(!aodExtraTracks)return iCount;
1905 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1906 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1907 if (!track) continue;
1909 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1910 if(!trackAOD)continue;
1911 Bool_t bGood = false;
1912 if(fFilterType == 0)bGood = true;
1913 else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1914 else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1915 if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1916 if(fRequireITSRefit==0){if((trackAOD->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
1917 if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1918 if(trackAOD->Pt()<fTrackPtCut) continue;
1919 if(fDebug) printf("pt extra track %.2f \n", trackAOD->Pt());
1920 list->Add(trackAOD);
1925 if(type==kTrackAODMCextra || type==kTrackAODMCextraonly) { //embed generator level particles
1926 AliAODEvent *aod = 0;
1927 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1928 else aod = AODEvent();
1932 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraMCparticles"));
1933 if(!aodExtraTracks)return iCount;
1934 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1935 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1936 AliAODMCParticle *partmc = dynamic_cast<AliAODMCParticle*> ((*aodExtraTracks)[it]);
1938 if(fDebug>10) printf("track %d does not exist\n",it);
1942 if(!partmc) continue;
1943 if(!partmc->IsPhysicalPrimary())continue;
1945 if (track->Pt()<fTrackPtCut) {
1946 if(fDebug>10) printf("track %d has too low pt %.2f\n",it,track->Pt());
1951 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*>((*aodExtraTracks)[it]);//(track);
1954 if(fDebug>10) printf("trackAOD %d does not exist\n",it);
1958 trackAOD->SetFlags(AliESDtrack::kEmbedded);
1959 trackAOD->SetFilterMap(fFilterMask);
1961 if(fDebug>10) printf("pt extra track %.2f \n", track->Pt());
1963 if(TMath::Abs(track->Eta())>fTrackEtaWindow) continue;
1964 if(track->Pt()<fTrackPtCut) continue;
1972 else if (type == kTrackKineAll||type == kTrackKineCharged){
1973 AliMCEvent* mcEvent = MCEvent();
1974 if(!mcEvent)return iCount;
1975 // we want to have alivpartilces so use get track
1976 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1977 if(!mcEvent->IsPhysicalPrimary(it))continue;
1978 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1979 if(type == kTrackKineAll){
1980 if(part->Pt()<fTrackPtCut)continue;
1984 else if(type == kTrackKineCharged){
1985 if(part->Particle()->GetPDG()->Charge()==0)continue;
1986 if(part->Pt()<fTrackPtCut)continue;
1992 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1993 AliAODEvent *aod = 0;
1994 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1995 else aod = AODEvent();
1996 if(!aod)return iCount;
1997 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1998 if(!tca)return iCount;
1999 for(int it = 0;it < tca->GetEntriesFast();++it){
2000 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
2001 if(!part->IsPhysicalPrimary())continue;
2002 if(type == kTrackAODMCAll){
2003 if(part->Pt()<fTrackPtCut)continue;
2007 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
2008 if(part->Charge()==0)continue;
2009 if(part->Pt()<fTrackPtCut)continue;
2010 if(kTrackAODMCCharged){
2014 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
2021 else if (type == kTrackAODMCHF){
2023 AliAODEvent *aod = 0;
2024 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
2025 else aod = AODEvent();
2026 if(!aod)return iCount;
2027 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
2028 if(!tca)return iCount;
2029 for(int it = 0;it < tca->GetEntriesFast();++it){
2030 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
2032 int partpdg= part->PdgCode();
2033 if(!part->IsPhysicalPrimary() && !IsBMeson(partpdg) && !IsDMeson(partpdg) )continue;
2035 if (IsBMeson(partpdg) || IsDMeson(partpdg)) {
2036 iCount+= AddDaughters( list , part,tca);
2040 if(part->Pt()<fTrackPtCut) continue;
2041 if(TMath::Abs(part->Eta())>fTrackEtaWindow) continue;
2042 if(part->Charge()==0) continue;
2044 if((part->Pt()>=fTrackPtCut) && (TMath::Abs(part->Eta())<=fTrackEtaWindow) && (part->Charge()!=0))list->Add(part);
2054 Int_t AliAnalysisTaskJetCluster::AddDaughters(TList * list, AliAODMCParticle *part, TClonesArray * tca){
2056 Int_t nDaugthers = part->GetNDaughters();
2057 for(Int_t i=0;i< nDaugthers;++i){
2058 AliAODMCParticle *partdaughter = (AliAODMCParticle*)(tca->At(i));
2059 if(!partdaughter) continue;
2060 if(partdaughter->Pt()<fTrackPtCut)continue;
2061 if(TMath::Abs(partdaughter->Eta())>fTrackEtaWindow)continue;
2062 if(partdaughter->Charge()==0)continue;
2064 if(!IsDMeson(partdaughter->PdgCode()) && !IsBMeson(partdaughter->PdgCode()) ){
2065 list->Add(partdaughter);
2068 else count+=AddDaughters(list,part,tca);
2075 void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() {
2078 AliInfo("Trying to connect to AliEn ...");
2079 TGrid::Connect("alien://");
2082 TFile *f = TFile::Open(fPathTrPtResolution.Data());
2086 TProfile *fProfPtPtSigma1PtGlobSt = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
2087 TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
2088 TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
2090 SetSmearResolution(kTRUE);
2091 SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD);
2096 void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() {
2099 AliInfo("Trying to connect to AliEn ...");
2100 TGrid::Connect("alien://");
2103 TFile *f = TFile::Open(fPathTrEfficiency.Data());
2106 TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
2107 TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
2108 TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
2110 SetDiceEfficiency(kTRUE);
2112 if(fChangeEfficiencyFraction>0.) {
2114 TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin");
2116 for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) {
2117 Double_t content = hEffPosGlobSt->GetBinContent(bin);
2118 hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction);
2121 SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
2125 SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
2129 void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
2132 // set mom res profiles
2135 if(fMomResH1) delete fMomResH1;
2136 if(fMomResH2) delete fMomResH2;
2137 if(fMomResH3) delete fMomResH3;
2139 fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
2140 fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
2141 fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
2145 void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
2147 // set tracking efficiency histos
2150 fhEffH1 = (TH1*)h1->Clone("fhEffH1");
2151 fhEffH2 = (TH1*)h2->Clone("fhEffH2");
2152 fhEffH3 = (TH1*)h3->Clone("fhEffH3");
2155 Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
2158 // Get smearing on generated momentum
2161 //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
2163 TProfile *fMomRes = 0x0;
2164 if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
2165 if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
2166 if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
2173 Double_t smear = 0.;
2176 if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
2177 if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
2178 if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
2182 Int_t bin = fMomRes->FindBin(pt);
2184 smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
2188 if(fMomRes) delete fMomRes;
2193 void AliAnalysisTaskJetCluster::FitMomentumResolution() {
2195 // Fit linear function on momentum resolution at high pT
2198 if(!fMomResH1Fit && fMomResH1) {
2199 fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
2200 fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
2201 fMomResH1Fit ->SetRange(5.,100.);
2204 if(!fMomResH2Fit && fMomResH2) {
2205 fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
2206 fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
2207 fMomResH2Fit ->SetRange(5.,100.);
2210 if(!fMomResH3Fit && fMomResH3) {
2211 fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
2212 fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
2213 fMomResH3Fit ->SetRange(5.,100.);
2219 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
2220 for(int i = 0; i < particles.GetEntries(); i++){
2221 AliVParticle *vp = (AliVParticle*)particles.At(i);
2222 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
2223 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
2224 jInp.set_user_index(i);
2225 inputParticles.push_back(jInp);
2234 bool AliAnalysisTaskJetCluster::IsBMeson(int pc){
2235 int bPdG[] = {511,521,10511,10521,513,523,10513,10523,20513,20523,20513,20523,515,525,531,
2236 10531,533,10533,20533,535,541,10541,543,10543,20543,545,551,10551,100551,
2237 110551,200551,210551,553,10553,20553,30553,100553,110553,120553,130553,200553,210553,220553,
2238 300553,9000533,9010553,555,10555,20555,100555,110555,120555,200555,557,100557};
2239 for(int i=0;i< (int)(sizeof(bPdG)/sizeof(int));++i) if(abs(pc) == bPdG[i]) return true;
2242 bool AliAnalysisTaskJetCluster::IsDMeson(int pc){
2243 int bPdG[] = {411,421,10411,10421,413,423,10413,10423,20431,20423,415,
2244 425,431,10431,433,10433,20433,435,441,10441,100441,443,10443,20443,
2245 100443,30443,9000443,9010443,9020443,445,100445};
2246 for(int i=0;i< (int)(sizeof(bPdG)/sizeof(int));++i) if(abs(pc) == bPdG[i]) return true;