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"
48 #include "AliAODHandler.h"
50 #include "AliAnalysisTaskJetResponseV2.h"
55 ClassImp(AliAnalysisTaskJetResponseV2)
57 AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2() :
64 fBackgroundBranch(""),
66 fOfflineTrgMask(AliVEvent::kAny),
79 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
80 fJetPtFractionMin(0.5),
88 fbJetsMismatch1(kTRUE),
89 fbJetsMismatch2(kTRUE),
95 fbJets3Branches(kFALSE),
96 fbJetsBeforeCut1(kTRUE),
97 fbJetsBeforeCut2(kTRUE),
98 fHistEvtSelection(0x0),
99 fHistJetSelection(0x0),
100 fh2JetSelection(0x0),
102 fhnJetsMismatch1(0x0),
103 fhnJetsMismatch2(0x0),
109 fhnJets3Branches(0x0),
110 fhnJetsBeforeCut1(0x0),
111 fhnJetsBeforeCut2(0x0)
113 // default Constructor
115 fJetBranchName[0] = "";
116 fJetBranchName[1] = "";
117 fJetBranchName[2] = "";
119 fListJets[0] = new TList;
120 fListJets[1] = new TList;
121 fListJets[2] = new TList;
124 AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2(const char *name) :
125 AliAnalysisTaskSE(name),
131 fBackgroundBranch(""),
133 fOfflineTrgMask(AliVEvent::kAny),
142 fNInputTracksMax(-1),
146 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
147 fJetPtFractionMin(0.5),
155 fbJetsMismatch1(kTRUE),
156 fbJetsMismatch2(kTRUE),
158 fbJetsDeltaPt(kTRUE),
162 fbJets3Branches(kFALSE),
163 fbJetsBeforeCut1(kTRUE),
164 fbJetsBeforeCut2(kTRUE),
165 fHistEvtSelection(0x0),
166 fHistJetSelection(0x0),
167 fh2JetSelection(0x0),
169 fhnJetsMismatch1(0x0),
170 fhnJetsMismatch2(0x0),
176 fhnJets3Branches(0x0),
177 fhnJetsBeforeCut1(0x0),
178 fhnJetsBeforeCut2(0x0)
182 fJetBranchName[0] = "";
183 fJetBranchName[1] = "";
184 fJetBranchName[2] = "";
186 fListJets[0] = new TList;
187 fListJets[1] = new TList;
188 fListJets[2] = new TList;
190 DefineOutput(1, TList::Class());
193 AliAnalysisTaskJetResponseV2::~AliAnalysisTaskJetResponseV2()
200 void AliAnalysisTaskJetResponseV2::SetBranchNames(const TString &branch1, const TString &branch2, const TString &branch3)
202 fJetBranchName[0] = branch1;
203 fJetBranchName[1] = branch2;
204 fJetBranchName[2] = branch3;
206 if(strlen(fJetBranchName[2].Data()) ) {
208 fbJets3Branches = kTRUE;
212 void AliAnalysisTaskJetResponseV2::Init()
215 // check for jet branches
216 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
217 AliError("Jet branch name not set.");
222 void AliAnalysisTaskJetResponseV2::UserCreateOutputObjects()
227 if(!fOutputList) fOutputList = new TList;
228 fOutputList->SetOwner(kTRUE);
230 Bool_t oldStatus = TH1::AddDirectoryStatus();
231 TH1::AddDirectory(kFALSE);
234 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
235 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
236 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
237 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
238 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
239 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
240 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
242 fHistJetSelection = new TH1I("fHistJetSelection", "jet selection", 8, -0.5, 7.5);
243 fHistJetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
244 fHistJetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
245 fHistJetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
246 fHistJetSelection->GetXaxis()->SetBinLabel(4,"not in list");
247 fHistJetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
248 fHistJetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
249 fHistJetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
250 fHistJetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
252 fh2JetSelection = new TH2F("fh2JetSelection", "jet selection", 8, -0.5, 7.5,100,0.,200.);
253 fh2JetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
254 fh2JetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
255 fh2JetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
256 fh2JetSelection->GetXaxis()->SetBinLabel(4,"not in list");
257 fh2JetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
258 fh2JetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
259 fh2JetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
260 fh2JetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
263 UInt_t entries = 0; // bit coded, see GetDimParams() below
264 UInt_t opt = 0; // bit coded, default (0) or high resolution (1)
267 entries = 1<<0 | 1<<1 | 1<<2 | 1<<26; // cent : nInpTrks : rp psi : pT hard bin
268 opt = 1<<0 | 1<<1; // centrality and nInpTrks in high resolution
269 fhnEvent = NewTHnSparseF("fhnEvent", entries, opt);
272 if(fbJetsMismatch1){ // real mismatch (no related rec jet found)
273 // cent : nInpTrks : rp bins : probe pt : probe eta : probe phi : probe area : pT hard bin
274 entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<8 | 1<<10 | 1<<12 | 1<<26;
275 opt = 1<<6 | 1<<8 | 1<<10;
276 fhnJetsMismatch1 = NewTHnSparseF("fhnJetsMismatch1", entries, opt);
279 if(fbJetsMismatch2){ // acceptance + fraction cut
280 // cent : nInpTrks : rp bins : jetPt(2x) : jetEta(2x) : deltaEta : deltaR : fraction : pT hard bin
281 entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<15 | 1<<17 | 1<<19 | 1<<26;
282 opt = 1<<6 | 1<<7 | 1<<8 | 1<<9;
283 fhnJetsMismatch2 = NewTHnSparseF("fhnJetsMismatch2", entries, opt);
289 // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : probe area: deltaPt : rp phi : rho : correction for RP : local rho
291 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;
293 // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : deltaPt : rp phi : rho : correction for RP | pT hard bin
294 entries = 1<<0 | 1<<1 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<14 | 1<<2 | 1<<22 | 1<<23 | 1<<26;
296 fhnJetsRp = NewTHnSparseF("fhnJetsRp", entries, opt);
299 // cent : nInpTrks : rp bins: deltaPt : jetPt(2x) : deltaArea : pT hard bin (hr delta pt)
301 entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<18 | 1<<26;
302 opt = 1<<1 | 1<<14 | 1<<6 | 1<<7 | 1<<18;
303 fhnJetsDeltaPt = NewTHnSparseF("fhnJetsDeltaPt", entries, opt);
306 // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) : pT hard bin (hr for eta)
308 entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9 | 1<<26;
309 opt = 1<<15 | 1<<8 | 1<<9;
310 fhnJetsEta = NewTHnSparseF("fhnJetsEta", entries, opt);
313 // cent : nInpTrks : rp bins : jetPt(2x) : jetPhi(2x) : deltaPt : deltaPhi : pT hard bin
315 entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<10 | 1<<11 | 1<<14 | 1<<16 | 1<<26;
317 fhnJetsPhi = NewTHnSparseF("fhnJetsPhi", entries, opt);
320 // 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)
322 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;
323 opt = 1<<18 | 1<<12 | 1<<13;
324 fhnJetsArea = NewTHnSparseF("fhnJetsArea", entries, opt);
327 // cent : nInpTrks : jetPt(3x) : deltaPt : delta : jetArea(3x) : fraction(2x) : deltaR(1x) : pT hard bin
329 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<<17 | 1<<26;
330 opt = 1<<6 | 1<<7 | 1<<27 | 1<<14 | 1<<28;
331 fhnJets3Branches = NewTHnSparseF("fhnJets3Branches", entries, opt);
339 // cent : nInpTrks : rp bins : fraction : jetPt(2x) : jetEta(2x) : jetPhi(2x) (low resolution) (with fraction, eta, phi, pt cuts possible)
340 if(fbJetsBeforeCut1){
341 entries = 1<<0 | 1<<1 | 1<<3 | 1<<19 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<10 | 1<<11;
343 fhnJetsBeforeCut1 = NewTHnSparseF("fhnJetsBeforeCut1", entries, opt);
346 // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) (low resolution)
347 if(fbJetsBeforeCut2){
348 entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9;
350 fhnJetsBeforeCut2 = NewTHnSparseF("fhnJetsBeforeCut2", entries, opt);
353 fOutputList->Add(fHistEvtSelection);
354 fOutputList->Add(fHistJetSelection);
355 fOutputList->Add(fh2JetSelection);
356 if(fbEvent) fOutputList->Add(fhnEvent);
357 if(fbJetsMismatch1) fOutputList->Add(fhnJetsMismatch1);
358 if(fbJetsMismatch2) fOutputList->Add(fhnJetsMismatch2);
359 if(fbJetsRp) fOutputList->Add(fhnJetsRp);
360 if(fbJetsDeltaPt) fOutputList->Add(fhnJetsDeltaPt);
361 if(fbJetsEta) fOutputList->Add(fhnJetsEta);
362 if(fbJetsPhi) fOutputList->Add(fhnJetsPhi);
363 if(fbJetsArea) fOutputList->Add(fhnJetsArea);
364 if(fbJets3Branches) fOutputList->Add(fhnJets3Branches);
365 if(fbJetsBeforeCut1) fOutputList->Add(fhnJetsBeforeCut1);
366 if(fbJetsBeforeCut2) fOutputList->Add(fhnJetsBeforeCut2);
368 // =========== Switch on Sumw2 for all histos ===========
369 for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
370 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
375 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
380 TH1::AddDirectory(oldStatus);
382 PostData(1, fOutputList);
385 void AliAnalysisTaskJetResponseV2::UserExec(Option_t *)
387 // load events, apply event cuts, then compare leading jets
389 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
390 AliError("Jet branch name not set.");
394 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
396 AliError("ESD not available");
398 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
399 // assume that the AOD is in the general output...
400 fAODOut = AODEvent();
402 fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
405 AliError("AOD not available");
409 if(fNonStdFile.Length()!=0){
410 // case that we have an AOD extension we need can fetch the jets from the extended output
411 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
412 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
414 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
418 // -- event selection --
419 fHistEvtSelection->Fill(1); // number of events before event selection
422 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
423 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
424 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
425 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
426 fHistEvtSelection->Fill(2);
427 PostData(1, fOutputList);
432 AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
433 Int_t nTracksPrim = primVtx->GetNContributors();
434 if ((nTracksPrim < fMinContribVtx) ||
435 (primVtx->GetZ() < fVtxZMin) ||
436 (primVtx->GetZ() > fVtxZMax) ){
437 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
438 fHistEvtSelection->Fill(3);
439 PostData(1, fOutputList);
443 // event class selection (from jet helper task)
444 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
445 if(fDebug) Printf("Event class %d", eventClass);
446 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
447 fHistEvtSelection->Fill(4);
448 PostData(1, fOutputList);
452 // centrality selection
453 AliCentrality *cent = 0x0;
454 Float_t centValue = 0.;
455 if(fESD) cent = fESD->GetCentrality();
456 if(cent) centValue = cent->GetCentralityPercentile("V0M");
457 if(fDebug) printf("centrality: %f\n", centValue);
458 if (centValue < fCentMin || centValue > fCentMax){
459 fHistEvtSelection->Fill(4);
460 PostData(1, fOutputList);
465 // multiplicity due to input tracks
466 Int_t nInputTracks = GetNInputTracks();
468 if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
469 fHistEvtSelection->Fill(5);
470 PostData(1, fOutputList);
475 fHistEvtSelection->Fill(0); // accepted events
476 // -- end event selection --
479 Double_t pthard = AliAnalysisTaskFastEmbedding::GetPtHard();
480 Int_t pthardbin = GetPtHardBin(pthard);
483 Double_t rp = AliAnalysisHelperJetTasks::ReactionPlane(kFALSE);
485 Double_t eventEntries[3] = { (Double_t)centValue, (Double_t)nInputTracks, rp };
486 fhnEvent->Fill(eventEntries);
491 AliAODJetEventBackground* externalBackground = 0;
492 if(!externalBackground&&fBackgroundBranch.Length()){
493 externalBackground = (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
494 if(!externalBackground && fAODOut) externalBackground = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
495 if(!externalBackground && fAODExtension) externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
496 //if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
499 if(externalBackground) rho = externalBackground->GetBackground(0);
503 TClonesArray *aodJets[3];
504 aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
505 if(!aodJets[0] && fAODOut) aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
506 if(!aodJets[0] && fAODExtension) aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
507 aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE version1
508 if(!aodJets[1] && fAODOut) aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet
509 if(!aodJets[1] && fAODExtension) aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));
510 if( strlen(fJetBranchName[2].Data()) ) {
511 aodJets[2] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[2].Data())); // in general: embedded jet + UE version2
512 if(!aodJets[2] && fAODOut) aodJets[2] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[2].Data())); // in general: embedded jet
513 if(!aodJets[2] && fAODExtension) aodJets[2] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[2].Data()));
514 if(aodJets[2]) fkNbranches=3;
515 if(fDebug>10) printf("3rd branch: %s\n",fJetBranchName[2].Data());
518 if(fDebug>10) printf("fkNbranches %d\n",fkNbranches);
519 for (Int_t iJetType = 0; iJetType < fkNbranches; iJetType++) {
520 fListJets[iJetType]->Clear();
521 if (!aodJets[iJetType]) {
522 if(fDebug) Printf("%s: no jet branch\n",fJetBranchName[iJetType].Data());
525 if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
527 for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
528 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
529 if (jet) fListJets[iJetType]->Add(jet);
531 fListJets[iJetType]->Sort();
536 static TArrayI aMatchIndex(fListJets[0]->GetEntries());
537 static TArrayF aPtFraction(fListJets[0]->GetEntries());
538 if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
539 if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
541 // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
542 AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
543 fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
544 aMatchIndex, aPtFraction, fDebug, fMatchMaxDist, fIsPbPb?1:2);
545 static TArrayI aMatchIndexv2(fListJets[0]->GetEntries());
546 static TArrayF aPtFractionv2(fListJets[0]->GetEntries());
547 if(aMatchIndexv2.GetSize()<fListJets[0]->GetEntries()) aMatchIndexv2.Set(fListJets[0]->GetEntries());
548 if(aPtFractionv2.GetSize()<fListJets[0]->GetEntries()) aPtFractionv2.Set(fListJets[0]->GetEntries());
549 if( strlen(fJetBranchName[2].Data()) ) {
550 AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
551 fListJets[2], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[2]->GetEntries()),
552 aMatchIndexv2, aPtFractionv2, fDebug, fMatchMaxDist, fIsPbPb?1:2);
555 // loop over matched jets
556 Int_t ir = -1; // index of matched reconstruced jet
557 Int_t ir2 = -1; // index of matched reconstruced jet
558 Float_t fraction = -1.;
559 Float_t fraction2 = -1.;
560 AliAODJet *jet[3] = { 0x0, 0x0, 0x0 };
561 Float_t jetEta[3] = { -990., -990., -990. };
562 Float_t jetPhi[3] = { -990., -990., -990. };
563 Float_t jetPt[3] = { -990., -990., -990. };
564 Float_t jetArea[3] = { -990., -990., -990. };
565 Float_t rpJet[3] = { -990., -990., -990. };
568 for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
569 ir = aMatchIndex[ig];
570 if(aMatchIndexv2.GetSize()>ig) ir2 = aMatchIndexv2[ig];
574 jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
575 if(ir>=0) jet[1] = (AliAODJet*)(fListJets[1]->At(ir));
577 if(ir2>=0) jet[2] = (AliAODJet*)(fListJets[2]->At(ir2));
580 for(Int_t i=0; i<fkNbranches; ++i){
587 if(i==1) rpBin = -990;
590 jetEta[i] = jet[i]->Eta();
591 jetPhi[i] = jet[i]->Phi();
592 jetPt[i] = GetPt(jet[i], i);
593 jetArea[i] = jet[i]->EffectiveAreaCharged();
594 rpJet[i] = TVector2::Phi_mpi_pi(rp-jetPhi[i]);
595 if(i==1) rpBin = AliAnalysisHelperJetTasks::GetPhiBin(TVector2::Phi_mpi_pi(rp-jetPhi[i]), 3);
598 fraction = aPtFraction[ig];
599 if(aPtFractionv2.GetSize()>ig) fraction2 = aPtFractionv2[ig];
603 fHistJetSelection->Fill(1); // all probe jets
604 if(jet[0]) fh2JetSelection->Fill(1.,jetPt[0]); // all probe jets
607 fHistJetSelection->Fill(2);
608 if(jet[0]) fh2JetSelection->Fill(2.,jetPt[0]);
611 if(!jet[0]) continue;
612 Double_t jetEntriesMismatch1[7] = {
613 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
614 (Double_t)jetPt[0], (Double_t)jetEta[0], (Double_t)jetPhi[0], (Double_t)jetArea[0]
616 fhnJetsMismatch1->Fill(jetEntriesMismatch1);
622 if(!jet[0] || !jet[1]){
623 fHistJetSelection->Fill(3);
624 if(jet[0]) fh2JetSelection->Fill(3.,jetPt[0]);
628 // look for distance to next rec jet
629 Float_t distNextJet = -0.01; // no neighbor
630 Float_t ptNextJet = -1.;
631 for(Int_t r=0; r<fListJets[1]->GetEntries(); ++r){
633 Float_t tmpDeltaR = jet[1]->DeltaR((AliAODJet*)fListJets[1]->At(r));
634 if(distNextJet<0. || distNextJet>tmpDeltaR){
635 distNextJet = tmpDeltaR;
636 if(fKeepJets) ptNextJet = ((AliAODJet*)fListJets[1]->At(r))->GetPtSubtracted(0);
637 else ptNextJet = ((AliAODJet*)fListJets[1]->At(r))->Pt();
643 //Float_t localRho = jetArea[1]>0. ? (jetPt[1]+rho*jetArea[1] - jetPt[0]) / jetArea[1] : 0.;
644 //Float_t relRho = rho>0. ? localRho / rho : 0.;
647 // calculate parameters of associated jets
648 /* from Leticia, kT clusterizer
649 Float_t par0[4] = { 0.00409, 0.01229, 0.05, 0.26 };
650 Float_t par1[4] = { -2.97035e-03, -2.03182e-03, -1.25702e-03, -9.95107e-04 };
651 Float_t par2[4] = { 1.02865e-01, 1.49039e-01, 1.53910e-01, 1.51109e-01 };
653 // own, from embedded tracks
654 Float_t par0[4] = { 0.02841, 0.05039, 0.09092, 0.24089 };
655 Float_t par1[4] = { -4.26725e-04, -1.15273e-03, -1.56827e-03, -3.08003e-03 };
656 Float_t par2[4] = { 4.95415e-02, 9.79538e-02, 1.32814e-01, 1.71743e-01 };
660 if(eventClass>0&&eventClass<4){
661 rpCorr = par0[eventClass-1] + par1[eventClass-1] * 2*TMath::Cos(rpJet[1]) + par2[eventClass-1] * 2*TMath::Cos(2*rpJet[1]);
662 rpCorr *= rho * jetArea[1];
665 Float_t delta = jetPt[2]-jetPt[1];
666 Float_t deltaPt = jetPt[1]-jetPt[0];
667 Float_t deltaEta = jetEta[1]-jetEta[0];
668 Float_t deltaPhi = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]);
669 Float_t deltaR = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
670 Float_t deltaArea = jetArea[1]-jetArea[0];
673 // fill thnsparse before acceptance cut
674 if(fbJetsBeforeCut1){
675 Double_t jetBeforeCutEntries1[10] = {
676 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
677 (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1],
678 (Double_t)jetPhi[0], (Double_t)jetPhi[1], (Double_t)fraction
680 fhnJetsBeforeCut1->Fill(jetBeforeCutEntries1);
683 if(fbJetsBeforeCut2){
684 Double_t jetBeforeCutEntries2[10] = {
685 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
686 (Double_t)jetPt[0], (Double_t)jetPt[1],
687 (Double_t)jetEta[0], (Double_t)jetEta[1],
688 (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR
690 fhnJetsBeforeCut2->Fill(jetBeforeCutEntries2);
694 Bool_t jetAccepted = kTRUE;
695 // minimum fraction required
696 if(fraction<fJetPtFractionMin){
697 fHistJetSelection->Fill(4);
698 fh2JetSelection->Fill(4.,jetPt[0]);
699 jetAccepted = kFALSE;
703 // jet acceptance + minimum pT check
704 if(jetEta[0]>fJetEtaMax || jetEta[0]<fJetEtaMin ||
705 jetEta[1]>fJetEtaMax || jetEta[1]<fJetEtaMin){
708 Printf("Jet not in eta acceptance.");
709 Printf("[0]: jet %d eta %.2f", ig, jetEta[0]);
710 Printf("[1]: jet %d eta %.2f", ir, jetEta[1]);
712 fHistJetSelection->Fill(5);
713 fh2JetSelection->Fill(5.,jetPt[0]);
714 jetAccepted = kFALSE;
718 if(jetPt[0] < fJetPtMin){
719 if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[0]);
720 fHistJetSelection->Fill(6);
721 fh2JetSelection->Fill(6.,jetPt[0]);
722 jetAccepted = kFALSE;
727 if(jet[1]->Trigger()&fJetTriggerExcludeMask){
728 fHistJetSelection->Fill(7);
729 fh2JetSelection->Fill(7.,jetPt[0]);
730 jetAccepted = kFALSE;
737 Double_t jetEntriesMismatch2[10] = {
738 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
739 (Double_t)jetPt[0], (Double_t)jetPt[1],
740 (Double_t)jetEta[0], (Double_t)jetEta[1],
741 (Double_t)deltaEta, (Double_t)deltaR,
744 fhnJetsMismatch2->Fill(jetEntriesMismatch2);
750 fHistJetSelection->Fill(0);
751 fh2JetSelection->Fill(0.,jetPt[0]);
755 Double_t jetEntriesRp[11] = {
756 (Double_t)centValue, (Double_t)nInputTracks, (Double_t) rp,
757 (Double_t)rpBin, (Double_t)rpJet[0], (Double_t)rpJet[1],
758 (Double_t)jetPt[0], (Double_t)deltaPt, (Double_t)rho, (Double_t)rpCorr,
761 fhnJetsRp->Fill(jetEntriesRp);
765 Double_t jetEntriesDeltaPt[8] = {
766 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
767 (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)deltaPt, (Double_t)deltaArea,
770 fhnJetsDeltaPt->Fill(jetEntriesDeltaPt);
774 Double_t jetEntriesEta[11] = {
775 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
776 (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1],
777 (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR,
780 fhnJetsEta->Fill(jetEntriesEta);
784 Double_t jetEntriesPhi[10] = {
785 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
786 (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetPhi[0], (Double_t)jetPhi[1],
787 (Double_t)deltaPt, (Double_t)deltaPhi,
790 fhnJetsPhi->Fill(jetEntriesPhi);
794 Double_t jetEntriesArea[14] = {
795 (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
796 (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetArea[0], (Double_t)jetArea[1],
797 (Double_t)deltaPt, (Double_t)deltaR, (Double_t)deltaArea,
798 (Double_t)fraction, (Double_t)distNextJet, (Double_t)ptNextJet,
801 fhnJetsArea->Fill(jetEntriesArea);
805 Double_t jetEntries3Branches[14] = {
806 (Double_t)centValue, (Double_t)nInputTracks,
807 (Double_t)jetPt[0], (Double_t)jetPt[1],
808 (Double_t)jetArea[0], (Double_t)jetArea[1],
809 (Double_t)deltaPt, (Double_t)fraction, (Double_t)pthardbin,
810 (Double_t)jetPt[2],(Double_t)delta,(Double_t)jetArea[2], (Double_t)fraction2, (Double_t)deltaR
812 fhnJets3Branches->Fill(jetEntries3Branches);
817 PostData(1, fOutputList);
820 void AliAnalysisTaskJetResponseV2::Terminate(const Option_t *)
822 // Draw result to the screen
823 // Called once at the end of the query
825 if (!GetOutputData(1))
829 Int_t AliAnalysisTaskJetResponseV2::GetNInputTracks()
831 // number of tracks in the event, obtained via jet finder
833 Int_t nInputTracks = 0;
835 TString jbname(fJetBranchName[1]);
836 //needs complete event, use jets without background subtraction
837 for(Int_t i=1; i<=3; ++i){
838 if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
841 if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
842 else if(jbname.Contains("AODMCextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
843 if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
844 else if(jbname.Contains("AODMCextra")) jbname.ReplaceAll("AODextra","AOD");
846 if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
847 TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
848 if(!tmpAODjets && fAODOut) tmpAODjets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(jbname.Data()));
849 if(!tmpAODjets && fAODExtension) tmpAODjets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(jbname.Data()));
851 Printf("Jet branch %s not found", jbname.Data());
852 Printf("AliAnalysisTaskJetResponseV2::GetNInputTracks FAILED");
856 for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
857 AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
859 TRefArray *trackList = jet->GetRefTracks();
860 Int_t nTracks = trackList->GetEntriesFast();
861 nInputTracks += nTracks;
862 if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
864 if(fDebug) Printf("---> input tracks: %d", nInputTracks);
869 THnSparse* AliAnalysisTaskJetResponseV2::NewTHnSparseF(const char* name, UInt_t entries, UInt_t opt)
871 // generate new THnSparseF, axes are defined in GetDimParams()
874 UInt_t tmp = entries;
877 tmp = tmp &~ -tmp; // clear lowest bit
880 TString hnTitle(name);
881 const Int_t dim = count;
888 while(c<dim && i<32){
890 Bool_t highres = opt&(1<<i);
892 GetDimParams(i, highres, label, nbins[c], xmin[c], xmax[c]);
893 hnTitle += Form(";%s",label.Data());
901 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
904 void AliAnalysisTaskJetResponseV2::GetDimParams(Int_t iEntry, Bool_t hr, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
906 // stores label and binning of axis for THnSparse
908 const Double_t pi = TMath::Pi();
913 label = "V0 centrality (%)";
927 label = "nb. of input tracks";
947 label = "event plane #psi";
961 label = "event plane bin";
970 if(iEntry==4)label = "#Delta#phi(RP-jet) (probe)";
971 if(iEntry==5)label = "#Delta#phi(RP-jet) (rec)";
981 if(iEntry==6)label = "probe p_{T} (GeV/c)";
982 if(iEntry==7 || iEntry==27)label = "rec p_{T} (GeV/c)";
997 if(iEntry==8)label = "probe #eta";
998 if(iEntry==9)label = "rec #eta";
1013 if(iEntry==10)label = "probe #phi";
1014 if(iEntry==11)label = "rec #phi";
1016 nbins = 90; // modulo 18 (sectors)
1030 if(iEntry==12)label = "probe area";
1031 if(iEntry==13 || iEntry==29)label = "rec area";
1045 if(iEntry==14) label = "#delta p_{T}";
1046 if(iEntry==28) label = "#delta";
1059 label = "#Delta#eta";
1073 label = "#Delta#phi";
1101 label = "#Deltaarea";
1130 label = "distance to closest rec jet";
1144 label = "p_{T} of closest rec jet";
1158 label = "abs. correction of #rho for RP";
1165 label = "local #rho";
1172 label = "local #rho / #rho";
1179 label = "p_{T,hard} bin";
1189 //____________________________________________________________________________
1190 Int_t AliAnalysisTaskJetResponseV2::GetPtHardBin(Double_t ptHard){
1192 // returns pt hard bin (for LHC10e14, LHC11a1x, LHC11a2x simulations)
1194 const Int_t nBins = 10;
1195 Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1198 while(bin<nBins-1 && binLimits[bin+1]<ptHard){
1205 //____________________________________________________________________________
1206 Double_t AliAnalysisTaskJetResponseV2::GetPt(AliAODJet *j, Int_t mode=0){
1208 // returns jet pt, also negative pt after background subtraction if available
1212 if(fKeepJets && mode>0){ // background subtracted pt, can be negative
1213 pt = j->GetPtSubtracted(0);
1216 pt = j->Pt(); // if negative, pt=0.01