1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 // task compares jets in two or three branches,
18 // written for analysis of jet embedding in HI events
29 #include "THnSparse.h"
34 #include "AliAnalysisTask.h"
35 #include "AliAnalysisManager.h"
37 #include "AliVEvent.h"
38 #include "AliESDEvent.h"
39 #include "AliESDInputHandler.h"
40 #include "AliCentrality.h"
41 #include "AliAnalysisHelperJetTasks.h"
42 #include "AliInputEventHandler.h"
43 #include "AliAODJetEventBackground.h"
44 #include "AliAnalysisTaskFastEmbedding.h"
46 #include "AliAODEvent.h"
47 #include "AliAODJet.h"
49 #include "AliAnalysisTaskJetResponseV2.h"
54 ClassImp(AliAnalysisTaskJetResponseV2)
56 AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2() :
60 fBackgroundBranch(""),
62 fOfflineTrgMask(AliVEvent::kAny),
75 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
76 fJetPtFractionMin(0.5),
84 fbJetsMismatch1(kTRUE),
85 fbJetsMismatch2(kTRUE),
91 fbJets3Branches(kFALSE),
92 fbJetsBeforeCut1(kTRUE),
93 fbJetsBeforeCut2(kTRUE),
94 fHistEvtSelection(0x0),
95 fHistJetSelection(0x0),
98 fhnJetsMismatch1(0x0),
99 fhnJetsMismatch2(0x0),
105 fhnJets3Branches(0x0),
106 fhnJetsBeforeCut1(0x0),
107 fhnJetsBeforeCut2(0x0)
109 // default Constructor
111 fJetBranchName[0] = "";
112 fJetBranchName[1] = "";
113 fJetBranchName[2] = "";
115 fListJets[0] = new TList;
116 fListJets[1] = new TList;
117 fListJets[2] = new TList;
120 AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2(const char *name) :
121 AliAnalysisTaskSE(name),
124 fBackgroundBranch(""),
126 fOfflineTrgMask(AliVEvent::kAny),
135 fNInputTracksMax(-1),
139 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
140 fJetPtFractionMin(0.5),
148 fbJetsMismatch1(kTRUE),
149 fbJetsMismatch2(kTRUE),
151 fbJetsDeltaPt(kTRUE),
155 fbJets3Branches(kFALSE),
156 fbJetsBeforeCut1(kTRUE),
157 fbJetsBeforeCut2(kTRUE),
158 fHistEvtSelection(0x0),
159 fHistJetSelection(0x0),
160 fh2JetSelection(0x0),
162 fhnJetsMismatch1(0x0),
163 fhnJetsMismatch2(0x0),
169 fhnJets3Branches(0x0),
170 fhnJetsBeforeCut1(0x0),
171 fhnJetsBeforeCut2(0x0)
175 fJetBranchName[0] = "";
176 fJetBranchName[1] = "";
177 fJetBranchName[2] = "";
179 fListJets[0] = new TList;
180 fListJets[1] = new TList;
181 fListJets[2] = new TList;
183 DefineOutput(1, TList::Class());
186 AliAnalysisTaskJetResponseV2::~AliAnalysisTaskJetResponseV2()
193 void AliAnalysisTaskJetResponseV2::SetBranchNames(const TString &branch1, const TString &branch2, const TString &branch3)
195 fJetBranchName[0] = branch1;
196 fJetBranchName[1] = branch2;
197 fJetBranchName[2] = branch3;
199 if(strlen(fJetBranchName[2].Data()) ) {
201 fbJets3Branches = kTRUE;
205 void AliAnalysisTaskJetResponseV2::Init()
208 // check for jet branches
209 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
210 AliError("Jet branch name not set.");
215 void AliAnalysisTaskJetResponseV2::UserCreateOutputObjects()
220 if(!fOutputList) fOutputList = new TList;
221 fOutputList->SetOwner(kTRUE);
223 Bool_t oldStatus = TH1::AddDirectoryStatus();
224 TH1::AddDirectory(kFALSE);
227 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
228 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
229 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
230 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
231 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
232 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
233 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
235 fHistJetSelection = new TH1I("fHistJetSelection", "jet selection", 8, -0.5, 7.5);
236 fHistJetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
237 fHistJetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
238 fHistJetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
239 fHistJetSelection->GetXaxis()->SetBinLabel(4,"not in list");
240 fHistJetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
241 fHistJetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
242 fHistJetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
243 fHistJetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
245 fh2JetSelection = new TH2F("fh2JetSelection", "jet selection", 8, -0.5, 7.5,100,0.,200.);
246 fh2JetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
247 fh2JetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
248 fh2JetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
249 fh2JetSelection->GetXaxis()->SetBinLabel(4,"not in list");
250 fh2JetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
251 fh2JetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
252 fh2JetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
253 fh2JetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
256 UInt_t entries = 0; // bit coded, see GetDimParams() below
257 UInt_t opt = 0; // bit coded, default (0) or high resolution (1)
260 entries = 1<<0 | 1<<1 | 1<<2 | 1<<26; // cent : nInpTrks : rp psi : pT hard bin
261 opt = 1<<0 | 1<<1; // centrality and nInpTrks in high resolution
262 fhnEvent = NewTHnSparseF("fhnEvent", entries, opt);
265 if(fbJetsMismatch1){ // real mismatch (no related rec jet found)
266 // cent : nInpTrks : rp bins : probe pt : probe eta : probe phi : probe area : pT hard bin
267 entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<8 | 1<<10 | 1<<12 | 1<<26;
268 opt = 1<<6 | 1<<8 | 1<<10;
269 fhnJetsMismatch1 = NewTHnSparseF("fhnJetsMismatch1", entries, opt);
272 if(fbJetsMismatch2){ // acceptance + fraction cut
273 // cent : nInpTrks : rp bins : jetPt(2x) : jetEta(2x) : deltaEta : deltaR : fraction : pT hard bin
274 entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<15 | 1<<17 | 1<<19 | 1<<26;
275 opt = 1<<6 | 1<<7 | 1<<8 | 1<<9;
276 fhnJetsMismatch2 = NewTHnSparseF("fhnJetsMismatch2", entries, opt);
282 // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : probe area: deltaPt : rp phi : rho : correction for RP : local rho
284 entries = 1<<0 | 1<<1 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<12 | 1<<14 | 1<<2 | 1<<22 | 1<<23 | 1<<24 | 1<<25;
286 // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : deltaPt : rp phi : rho : correction for RP | pT hard bin
287 entries = 1<<0 | 1<<1 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<14 | 1<<2 | 1<<22 | 1<<23 | 1<<26;
289 fhnJetsRp = NewTHnSparseF("fhnJetsRp", entries, opt);
292 // cent : nInpTrks : rp bins: deltaPt : jetPt(2x) : deltaArea : pT hard bin (hr delta pt)
294 entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<18 | 1<<26;
295 opt = 1<<1 | 1<<14 | 1<<6 | 1<<7 | 1<<18;
296 fhnJetsDeltaPt = NewTHnSparseF("fhnJetsDeltaPt", entries, opt);
299 // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) : pT hard bin (hr for eta)
301 entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9 | 1<<26;
302 opt = 1<<15 | 1<<8 | 1<<9;
303 fhnJetsEta = NewTHnSparseF("fhnJetsEta", entries, opt);
306 // cent : nInpTrks : rp bins : jetPt(2x) : jetPhi(2x) : deltaPt : deltaPhi : pT hard bin
308 entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<10 | 1<<11 | 1<<14 | 1<<16 | 1<<26;
310 fhnJetsPhi = NewTHnSparseF("fhnJetsPhi", entries, opt);
313 // cent : nInpTrks : rp bins : deltaArea : jetArea(2x) : deltaR : fraction : distance next rec jet : pT next jet : deltaPt : jetPt(2x) : pT hard bin (hr for area)
315 entries = 1<<0 | 1<<1 | 1<<3 | 1<<18 | 1<<12 | 1<<13 | 1<<17 | 1<<19 | 1<<20 | 1<<21 | 1<<14 | 1<<6 | 1<<7 | 1<<26;
316 opt = 1<<18 | 1<<12 | 1<<13;
317 fhnJetsArea = NewTHnSparseF("fhnJetsArea", entries, opt);
320 // cent : nInpTrks : jetPt(3x) : deltaPt : delta : jetArea(3x) : fraction(2x) : pT hard bin
322 entries = 1<<0 | 1<<1 | 1<<6 | 1<<7 | 1<<27 | 1<<14 | 1<<28 | 1<<12 | 1<<13 | 1<<29 | 1<<19 | 1<<30 | 1<<26;
323 opt = 1<<6 | 1<<7 | 1<<27 | 1<<14 | 1<<28;
324 fhnJets3Branches = NewTHnSparseF("fhnJets3Branches", entries, opt);
332 // cent : nInpTrks : rp bins : fraction : jetPt(2x) : jetEta(2x) : jetPhi(2x) (low resolution) (with fraction, eta, phi, pt cuts possible)
333 if(fbJetsBeforeCut1){
334 entries = 1<<0 | 1<<1 | 1<<3 | 1<<19 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<10 | 1<<11;
336 fhnJetsBeforeCut1 = NewTHnSparseF("fhnJetsBeforeCut1", entries, opt);
339 // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) (low resolution)
340 if(fbJetsBeforeCut2){
341 entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9;
343 fhnJetsBeforeCut2 = NewTHnSparseF("fhnJetsBeforeCut2", entries, opt);
346 fOutputList->Add(fHistEvtSelection);
347 fOutputList->Add(fHistJetSelection);
348 fOutputList->Add(fh2JetSelection);
349 if(fbEvent) fOutputList->Add(fhnEvent);
350 if(fbJetsMismatch1) fOutputList->Add(fhnJetsMismatch1);
351 if(fbJetsMismatch2) fOutputList->Add(fhnJetsMismatch2);
352 if(fbJetsRp) fOutputList->Add(fhnJetsRp);
353 if(fbJetsDeltaPt) fOutputList->Add(fhnJetsDeltaPt);
354 if(fbJetsEta) fOutputList->Add(fhnJetsEta);
355 if(fbJetsPhi) fOutputList->Add(fhnJetsPhi);
356 if(fbJetsArea) fOutputList->Add(fhnJetsArea);
357 if(fbJets3Branches) fOutputList->Add(fhnJets3Branches);
358 if(fbJetsBeforeCut1) fOutputList->Add(fhnJetsBeforeCut1);
359 if(fbJetsBeforeCut2) fOutputList->Add(fhnJetsBeforeCut2);
361 // =========== Switch on Sumw2 for all histos ===========
362 for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
363 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
368 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
373 TH1::AddDirectory(oldStatus);
375 PostData(1, fOutputList);
378 void AliAnalysisTaskJetResponseV2::UserExec(Option_t *)
380 // load events, apply event cuts, then compare leading jets
382 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
383 AliError("Jet branch name not set.");
387 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
389 AliError("ESD not available");
390 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
392 fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
395 AliError("AOD not available");
399 // -- event selection --
400 fHistEvtSelection->Fill(1); // number of events before event selection
403 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
404 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
405 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
406 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
407 fHistEvtSelection->Fill(2);
408 PostData(1, fOutputList);
413 AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
414 Int_t nTracksPrim = primVtx->GetNContributors();
415 if ((nTracksPrim < fMinContribVtx) ||
416 (primVtx->GetZ() < fVtxZMin) ||
417 (primVtx->GetZ() > fVtxZMax) ){
418 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
419 fHistEvtSelection->Fill(3);
420 PostData(1, fOutputList);
424 // event class selection (from jet helper task)
425 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
426 if(fDebug) Printf("Event class %d", eventClass);
427 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
428 fHistEvtSelection->Fill(4);
429 PostData(1, fOutputList);
433 // centrality selection
434 AliCentrality *cent = 0x0;
435 Float_t centValue = 0.;
436 if(fESD) cent = fESD->GetCentrality();
437 if(cent) centValue = cent->GetCentralityPercentile("V0M");
438 if(fDebug) printf("centrality: %f\n", centValue);
439 if (centValue < fCentMin || centValue > fCentMax){
440 fHistEvtSelection->Fill(4);
441 PostData(1, fOutputList);
446 // multiplicity due to input tracks
447 Int_t nInputTracks = GetNInputTracks();
449 if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
450 fHistEvtSelection->Fill(5);
451 PostData(1, fOutputList);
456 fHistEvtSelection->Fill(0); // accepted events
457 // -- end event selection --
460 Double_t pthard = AliAnalysisTaskFastEmbedding::GetPtHard();
461 Int_t pthardbin = GetPtHardBin(pthard);
464 Double_t rp = AliAnalysisHelperJetTasks::ReactionPlane(kFALSE);
466 Double_t eventEntries[3] = { (Double_t)centValue, (Double_t)nInputTracks, rp };
467 fhnEvent->Fill(eventEntries);
472 AliAODJetEventBackground* externalBackground = 0;
473 if(!externalBackground&&fBackgroundBranch.Length()){
474 externalBackground = (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
475 //if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
478 if(externalBackground)rho = externalBackground->GetBackground(0);
482 TClonesArray *aodJets[3];
483 aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
484 aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE version1
485 if( strlen(fJetBranchName[2].Data()) ) {
486 aodJets[2] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[2].Data())); // in general: embedded jet + UE version2
490 for (Int_t iJetType = 0; iJetType < fkNbranches; iJetType++) {
491 fListJets[iJetType]->Clear();
492 if (!aodJets[iJetType]) continue;
494 if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
496 for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
497 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
498 if (jet) fListJets[iJetType]->Add(jet);
500 fListJets[iJetType]->Sort();
505 static TArrayI aMatchIndex(fListJets[0]->GetEntries());
506 static TArrayF aPtFraction(fListJets[0]->GetEntries());
507 if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
508 if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
510 // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
511 AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
512 fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
513 aMatchIndex, aPtFraction, fDebug, fMatchMaxDist, fIsPbPb?1:2);
514 static TArrayI aMatchIndexv2(fListJets[0]->GetEntries());
515 static TArrayF aPtFractionv2(fListJets[0]->GetEntries());
516 if(aMatchIndexv2.GetSize()<fListJets[0]->GetEntries()) aMatchIndexv2.Set(fListJets[0]->GetEntries());
517 if(aPtFractionv2.GetSize()<fListJets[0]->GetEntries()) aPtFractionv2.Set(fListJets[0]->GetEntries());
518 if( strlen(fJetBranchName[2].Data()) ) {
519 AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
520 fListJets[2], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[2]->GetEntries()),
521 aMatchIndexv2, aPtFractionv2, fDebug, fMatchMaxDist, fIsPbPb?1:2);
524 // loop over matched jets
525 Int_t ir = -1; // index of matched reconstruced jet
526 Int_t ir2 = -1; // index of matched reconstruced jet
527 Float_t fraction = -1.;
528 Float_t fraction2 = -1.;
529 AliAODJet *jet[3] = { 0x0, 0x0, 0x0 };
530 Float_t jetEta[3] = { -990., -990., -990. };
531 Float_t jetPhi[3] = { -990., -990., -990. };
532 Float_t jetPt[3] = { -990., -990., -990. };
533 Float_t jetArea[3] = { -990., -990., -990. };
534 Float_t rpJet[3] = { -990., -990., -990. };
537 for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
538 ir = aMatchIndex[ig];
539 if(aMatchIndexv2.GetSize()>ig) ir2 = aMatchIndexv2[ig];
543 jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
544 if(ir>=0) jet[1] = (AliAODJet*)(fListJets[1]->At(ir));
546 if(ir2>=0) jet[2] = (AliAODJet*)(fListJets[2]->At(ir2));
549 for(Int_t i=0; i<fkNbranches; ++i){
556 if(i==1) rpBin = -990;
559 jetEta[i] = jet[i]->Eta();
560 jetPhi[i] = jet[i]->Phi();
561 jetPt[i] = GetPt(jet[i], i);
562 jetArea[i] = jet[i]->EffectiveAreaCharged();
563 rpJet[i] = TVector2::Phi_mpi_pi(rp-jetPhi[i]);
564 if(i==1) rpBin = AliAnalysisHelperJetTasks::GetPhiBin(TVector2::Phi_mpi_pi(rp-jetPhi[i]), 3);
567 fraction = aPtFraction[ig];
568 if(aPtFractionv2.GetSize()>ig) fraction2 = aPtFractionv2[ig];
572 fHistJetSelection->Fill(1); // all probe jets
573 if(jet[0]) fh2JetSelection->Fill(1.,jetPt[0]); // all probe jets
576 fHistJetSelection->Fill(2);
577 if(jet[0]) fh2JetSelection->Fill(2.,jetPt[0]);
580 if(!jet[0]) continue;
581 Double_t jetEntriesMismatch1[7] = {
582 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
583 (Double_t)jetPt[0], (Double_t)jetEta[0], (Double_t)jetPhi[0], (Double_t)jetArea[0]
585 fhnJetsMismatch1->Fill(jetEntriesMismatch1);
591 if(!jet[0] || !jet[1]){
592 fHistJetSelection->Fill(3);
593 if(jet[0]) fh2JetSelection->Fill(3.,jetPt[0]);
597 // look for distance to next rec jet
598 Float_t distNextJet = -0.01; // no neighbor
599 Float_t ptNextJet = -1.;
600 for(Int_t r=0; r<fListJets[1]->GetEntries(); ++r){
602 Float_t tmpDeltaR = jet[1]->DeltaR((AliAODJet*)fListJets[1]->At(r));
603 if(distNextJet<0. || distNextJet>tmpDeltaR){
604 distNextJet = tmpDeltaR;
605 if(fKeepJets) ptNextJet = ((AliAODJet*)fListJets[1]->At(r))->GetPtSubtracted(0);
606 else ptNextJet = ((AliAODJet*)fListJets[1]->At(r))->Pt();
612 //Float_t localRho = jetArea[1]>0. ? (jetPt[1]+rho*jetArea[1] - jetPt[0]) / jetArea[1] : 0.;
613 //Float_t relRho = rho>0. ? localRho / rho : 0.;
616 // calculate parameters of associated jets
617 /* from Leticia, kT clusterizer
618 Float_t par0[4] = { 0.00409, 0.01229, 0.05, 0.26 };
619 Float_t par1[4] = { -2.97035e-03, -2.03182e-03, -1.25702e-03, -9.95107e-04 };
620 Float_t par2[4] = { 1.02865e-01, 1.49039e-01, 1.53910e-01, 1.51109e-01 };
622 // own, from embedded tracks
623 Float_t par0[4] = { 0.02841, 0.05039, 0.09092, 0.24089 };
624 Float_t par1[4] = { -4.26725e-04, -1.15273e-03, -1.56827e-03, -3.08003e-03 };
625 Float_t par2[4] = { 4.95415e-02, 9.79538e-02, 1.32814e-01, 1.71743e-01 };
629 if(eventClass>0&&eventClass<4){
630 rpCorr = par0[eventClass-1] + par1[eventClass-1] * 2*TMath::Cos(rpJet[1]) + par2[eventClass-1] * 2*TMath::Cos(2*rpJet[1]);
631 rpCorr *= rho * jetArea[1];
634 Float_t delta = jetPt[2]-jetPt[1];
635 Float_t deltaPt = jetPt[1]-jetPt[0];
636 Float_t deltaEta = jetEta[1]-jetEta[0];
637 Float_t deltaPhi = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]);
638 Float_t deltaR = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
639 Float_t deltaArea = jetArea[1]-jetArea[0];
642 // fill thnsparse before acceptance cut
643 if(fbJetsBeforeCut1){
644 Double_t jetBeforeCutEntries1[10] = {
645 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
646 (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1],
647 (Double_t)jetPhi[0], (Double_t)jetPhi[1], (Double_t)fraction
649 fhnJetsBeforeCut1->Fill(jetBeforeCutEntries1);
652 if(fbJetsBeforeCut2){
653 Double_t jetBeforeCutEntries2[10] = {
654 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
655 (Double_t)jetPt[0], (Double_t)jetPt[1],
656 (Double_t)jetEta[0], (Double_t)jetEta[1],
657 (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR
659 fhnJetsBeforeCut2->Fill(jetBeforeCutEntries2);
663 Bool_t jetAccepted = kTRUE;
664 // minimum fraction required
665 if(fraction<fJetPtFractionMin){
666 fHistJetSelection->Fill(4);
667 fh2JetSelection->Fill(4.,jetPt[0]);
668 jetAccepted = kFALSE;
672 // jet acceptance + minimum pT check
673 if(jetEta[0]>fJetEtaMax || jetEta[0]<fJetEtaMin ||
674 jetEta[1]>fJetEtaMax || jetEta[1]<fJetEtaMin){
677 Printf("Jet not in eta acceptance.");
678 Printf("[0]: jet %d eta %.2f", ig, jetEta[0]);
679 Printf("[1]: jet %d eta %.2f", ir, jetEta[1]);
681 fHistJetSelection->Fill(5);
682 fh2JetSelection->Fill(5.,jetPt[0]);
683 jetAccepted = kFALSE;
687 if(jetPt[0] < fJetPtMin){
688 if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[0]);
689 fHistJetSelection->Fill(6);
690 fh2JetSelection->Fill(6.,jetPt[0]);
691 jetAccepted = kFALSE;
696 if(jet[1]->Trigger()&fJetTriggerExcludeMask){
697 fHistJetSelection->Fill(7);
698 fh2JetSelection->Fill(7.,jetPt[0]);
699 jetAccepted = kFALSE;
706 Double_t jetEntriesMismatch2[10] = {
707 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
708 (Double_t)jetPt[0], (Double_t)jetPt[1],
709 (Double_t)jetEta[0], (Double_t)jetEta[1],
710 (Double_t)deltaEta, (Double_t)deltaR,
713 fhnJetsMismatch2->Fill(jetEntriesMismatch2);
719 fHistJetSelection->Fill(0);
720 fh2JetSelection->Fill(0.,jetPt[0]);
724 Double_t jetEntriesRp[11] = {
725 (Double_t)centValue, (Double_t)nInputTracks, (Double_t) rp,
726 (Double_t)rpBin, (Double_t)rpJet[0], (Double_t)rpJet[1],
727 (Double_t)jetPt[0], (Double_t)deltaPt, (Double_t)rho, (Double_t)rpCorr,
730 fhnJetsRp->Fill(jetEntriesRp);
734 Double_t jetEntriesDeltaPt[8] = {
735 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
736 (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)deltaPt, (Double_t)deltaArea,
739 fhnJetsDeltaPt->Fill(jetEntriesDeltaPt);
743 Double_t jetEntriesEta[11] = {
744 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
745 (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1],
746 (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR,
749 fhnJetsEta->Fill(jetEntriesEta);
753 Double_t jetEntriesPhi[10] = {
754 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
755 (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetPhi[0], (Double_t)jetPhi[1],
756 (Double_t)deltaPt, (Double_t)deltaPhi,
759 fhnJetsPhi->Fill(jetEntriesPhi);
763 Double_t jetEntriesArea[14] = {
764 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
765 (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetArea[0], (Double_t)jetArea[1],
766 (Double_t)deltaPt, (Double_t)deltaR, (Double_t)deltaArea,
767 (Double_t)fraction, (Double_t)distNextJet, (Double_t)ptNextJet,
770 fhnJetsArea->Fill(jetEntriesArea);
774 Double_t jetEntries3Branches[13] = {
775 (Double_t)centValue, (Double_t)nInputTracks,
776 (Double_t)jetPt[0], (Double_t)jetPt[1],
777 (Double_t)jetArea[0], (Double_t)jetArea[1],
778 (Double_t)deltaPt, (Double_t)fraction, (Double_t)pthardbin,
779 (Double_t)jetPt[2],(Double_t)delta,(Double_t)jetArea[2], (Double_t)fraction2
781 fhnJets3Branches->Fill(jetEntries3Branches);
786 PostData(1, fOutputList);
789 void AliAnalysisTaskJetResponseV2::Terminate(const Option_t *)
791 // Draw result to the screen
792 // Called once at the end of the query
794 if (!GetOutputData(1))
798 Int_t AliAnalysisTaskJetResponseV2::GetNInputTracks()
800 // number of tracks in the event, obtained via jet finder
802 Int_t nInputTracks = 0;
804 TString jbname(fJetBranchName[1]);
805 //needs complete event, use jets without background subtraction
806 for(Int_t i=1; i<=3; ++i){
807 if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
810 if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
811 if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
813 if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
814 TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
816 Printf("Jet branch %s not found", jbname.Data());
817 Printf("AliAnalysisTaskJetResponseV2::GetNInputTracks FAILED");
821 for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
822 AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
824 TRefArray *trackList = jet->GetRefTracks();
825 Int_t nTracks = trackList->GetEntriesFast();
826 nInputTracks += nTracks;
827 if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
829 if(fDebug) Printf("---> input tracks: %d", nInputTracks);
834 THnSparse* AliAnalysisTaskJetResponseV2::NewTHnSparseF(const char* name, UInt_t entries, UInt_t opt)
836 // generate new THnSparseF, axes are defined in GetDimParams()
839 UInt_t tmp = entries;
842 tmp = tmp &~ -tmp; // clear lowest bit
845 TString hnTitle(name);
846 const Int_t dim = count;
853 while(c<dim && i<32){
855 Bool_t highres = opt&(1<<i);
857 GetDimParams(i, highres, label, nbins[c], xmin[c], xmax[c]);
858 hnTitle += Form(";%s",label.Data());
866 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
869 void AliAnalysisTaskJetResponseV2::GetDimParams(Int_t iEntry, Bool_t hr, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
871 // stores label and binning of axis for THnSparse
873 const Double_t pi = TMath::Pi();
878 label = "V0 centrality (%)";
892 label = "nb. of input tracks";
912 label = "event plane #psi";
926 label = "event plane bin";
935 if(iEntry==4)label = "#Delta#phi(RP-jet) (probe)";
936 if(iEntry==5)label = "#Delta#phi(RP-jet) (rec)";
946 if(iEntry==6)label = "probe p_{T} (GeV/c)";
947 if(iEntry==7 || iEntry==27)label = "rec p_{T} (GeV/c)";
962 if(iEntry==8)label = "probe #eta";
963 if(iEntry==9)label = "rec #eta";
978 if(iEntry==10)label = "probe #phi";
979 if(iEntry==11)label = "rec #phi";
981 nbins = 90; // modulo 18 (sectors)
995 if(iEntry==12)label = "probe area";
996 if(iEntry==13 || iEntry==29)label = "rec area";
1010 if(iEntry==14) label = "#delta p_{T}";
1011 if(iEntry==28) label = "#delta";
1024 label = "#Delta#eta";
1038 label = "#Delta#phi";
1066 label = "#Deltaarea";
1095 label = "distance to closest rec jet";
1109 label = "p_{T} of closest rec jet";
1123 label = "abs. correction of #rho for RP";
1130 label = "local #rho";
1137 label = "local #rho / #rho";
1144 label = "p_{T,hard} bin";
1154 //____________________________________________________________________________
1155 Int_t AliAnalysisTaskJetResponseV2::GetPtHardBin(Double_t ptHard){
1157 // returns pt hard bin (for LHC10e14, LHC11a1x, LHC11a2x simulations)
1159 const Int_t nBins = 10;
1160 Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1163 while(bin<nBins-1 && binLimits[bin+1]<ptHard){
1170 //____________________________________________________________________________
1171 Double_t AliAnalysisTaskJetResponseV2::GetPt(AliAODJet *j, Int_t mode=0){
1173 // returns jet pt, also negative pt after background subtraction if available
1177 if(fKeepJets && mode>0){ // background subtracted pt, can be negative
1178 pt = j->GetPtSubtracted(0);
1181 pt = j->Pt(); // if negative, pt=0.01