]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/AliAnalysisTaskJetBackgroundSubtract.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add option to use pTdet instead of pTgen as denominator of (pTgen - pTdet)/pT
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskJetBackgroundSubtract.cxx
CommitLineData
fb10c4b8 1// **************************************
2// Task used for the correction of determiantion of reconstructed jet spectra
3// Compares input (gen) and output (rec) jets
4// *******************************************
5
6
7/**************************************************************************
8 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9 * *
10 * Author: The ALICE Off-line Project. *
11 * Contributors are mentioned in the code where appropriate. *
12 * *
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 **************************************************************************/
21
22
23#include <TROOT.h>
24#include <TH1F.h>
d9c0c425 25#include <TH2F.h>
fb10c4b8 26#include <THnSparse.h>
27#include <TSystem.h>
f0659f11 28#include <TObjArray.h>
fb10c4b8 29#include <TInterpreter.h>
30#include <TList.h>
31#include <TLorentzVector.h>
d9c0c425 32#include <TRefArray.h>
fb10c4b8 33#include "TDatabasePDG.h"
34
35#include "AliAnalysisTaskJetBackgroundSubtract.h"
36#include "AliAnalysisManager.h"
37#include "AliAODHandler.h"
38#include "AliAODTrack.h"
39#include "AliAODJet.h"
40#include "AliAODEvent.h"
41#include "AliInputEventHandler.h"
42#include "AliAODJetEventBackground.h"
43
44
45ClassImp(AliAnalysisTaskJetBackgroundSubtract)
46
47AliAnalysisTaskJetBackgroundSubtract::~AliAnalysisTaskJetBackgroundSubtract(){
f0659f11 48 //
49 // destructor
50 //
fb10c4b8 51 delete fJBArray;
52 delete fOutJetArrayList;
53 delete fInJetArrayList;
54}
55
56AliAnalysisTaskJetBackgroundSubtract::AliAnalysisTaskJetBackgroundSubtract():
57 AliAnalysisTaskSE(),
58 fAODOut(0x0),
59 fAODIn(0x0),
60 fAODExtension(0x0),
f0659f11 61 fJBArray(0),
fb10c4b8 62 fBackgroundBranch(""),
63 fNonStdFile(""),
64 fReplaceString1("B0"),
65 fReplaceString2("B%d"),
d73e1768 66 fSubtraction(k4Area),
bea24d2d 67 fKeepJets(kFALSE),
e77ab2dc 68 fExclude2Leading(kTRUE),
fb10c4b8 69 fInJetArrayList(0x0),
70 fOutJetArrayList(0x0),
d73e1768 71 fh2CentvsRho(0x0),
72 fh2CentvsSigma(0x0),
82ac956f 73 fh2MultvsRho(0x0),
74 fh2MultvsSigma(0x0),
d73e1768 75 fh2ShiftEta(0x0),
76 fh2ShiftPhi(0x0),
77 fh2ShiftEtaLeading(0x0),
78 fh2ShiftPhiLeading(0x0),
fb10c4b8 79 fHistList(0x0)
80{
81
82}
83
84AliAnalysisTaskJetBackgroundSubtract::AliAnalysisTaskJetBackgroundSubtract(const char* name):
85
86 AliAnalysisTaskSE(name),
87 fAODOut(0x0),
88 fAODIn(0x0),
89 fAODExtension(0x0),
f0659f11 90 fJBArray(0),
fb10c4b8 91 fBackgroundBranch(""),
92 fNonStdFile(""),
93 fReplaceString1("B0"),
94 fReplaceString2("B%d"),
d73e1768 95 fSubtraction(k4Area),
bea24d2d 96 fKeepJets(kFALSE),
e77ab2dc 97 fExclude2Leading(kTRUE),
fb10c4b8 98 fInJetArrayList(0x0),
99 fOutJetArrayList(0x0),
d73e1768 100 fh2CentvsRho(0x0),
101 fh2CentvsSigma(0x0),
82ac956f 102 fh2MultvsRho(0x0),
103 fh2MultvsSigma(0x0),
d73e1768 104 fh2ShiftEta(0x0),
105 fh2ShiftPhi(0x0),
106 fh2ShiftEtaLeading(0x0),
107 fh2ShiftPhiLeading(0x0),
fb10c4b8 108 fHistList(0x0)
109{
110 DefineOutput(1, TList::Class());
111}
112
113
114
115Bool_t AliAnalysisTaskJetBackgroundSubtract::Notify()
116{
f0659f11 117 //
118 // exec with every new file
fb10c4b8 119 //
1c21a0e7 120 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
121
122 ResetOutJets();
123
124 // Now we also have the Input Event Available! Fillvthe pointers in the list
125
126 fInJetArrayList->Clear();
127 fOutJetArrayList->Clear();
128
f0659f11 129 for(int iJB = 0;iJB<(fJBArray?fJBArray->GetEntries():0);iJB++){
1c21a0e7 130 TObjString *ostr = (TObjString*)fJBArray->At(iJB);
a18eedbb 131
132
1c21a0e7 133 TClonesArray* jarray = 0;
134 if(!jarray&&fAODOut){
135 jarray = (TClonesArray*)(fAODOut->FindListObject(ostr->GetString().Data()));
136 }
137 if(!jarray&&fAODExtension){
138 jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(ostr->GetString().Data()));
139 }
140 if(!jarray&&fAODIn){
141 jarray = (TClonesArray*)(fAODIn->FindListObject(ostr->GetString().Data()));
142 }
143
144 if(!jarray){
145 if(fDebug){
146 Printf("%s:%d Input jet branch %s not found",(char*)__FILE__,__LINE__,ostr->GetString().Data());
147 }
148 continue;
149 }
150
151 TString newName(ostr->GetString().Data());
152 newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
153 TClonesArray* jarrayOut = 0;
ea950747 154 if(newName.CompareTo(ostr->GetString())==0){
155 Printf("%s:%d Input and output branch would have the same name, skipping %s ",(char*)__FILE__,__LINE__,ostr->GetString().Data());
156 continue;
157 }
158
1c21a0e7 159 if(!jarrayOut&&fAODOut){
160 jarrayOut = (TClonesArray*)(fAODOut->FindListObject(newName.Data()));
161 }
162 if(!jarrayOut&&fAODExtension){
163 jarrayOut = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(newName.Data()));
164 }
165
166 if(!jarrayOut){
167 if(fDebug){
168 Printf("%s:%d Output jet branch %s not found",(char*)__FILE__,__LINE__,newName.Data());
169 PrintAODContents();
170 }
171 continue;
172 }
173 if(jarrayOut&&jarray){
174 fOutJetArrayList->Add(jarrayOut);
175 fInJetArrayList->Add(jarray);
176 }
177 }
fb10c4b8 178 return kTRUE;
179}
180
181void AliAnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects()
182{
183
184 //
185 // Create the output container
186 //
187 // Connect the AOD
188
fb10c4b8 189 if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects() \n");
750e22e8 190
191
192
fb10c4b8 193 if(fNonStdFile.Length()!=0){
a18eedbb 194
fb10c4b8 195 // case that we have an AOD extension we need to fetch the jets from the extended output
196 // we identifay the extension aod event by looking for the branchname
197 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
8ddcf605 198
199 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
fb10c4b8 200
201 if(!fAODExtension){
202 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
203 }
204 }
205 fAODOut = AODEvent();
fb10c4b8 206
1c21a0e7 207 // usually we do not have the input already here
fb10c4b8 208
209 if(!fInJetArrayList)fInJetArrayList =new TList();
210 if(!fOutJetArrayList)fOutJetArrayList =new TList();
211
f0659f11 212 for(int iJB = 0;iJB<(fJBArray?fJBArray->GetEntries():0);iJB++){
fb10c4b8 213 TObjString *ostr = (TObjString*)fJBArray->At(iJB);
1c21a0e7 214 TString newName(ostr->GetString().Data());
fb10c4b8 215 if(!newName.Contains(fReplaceString1.Data())){
216 Printf("%s:%d cannot replace string %s in %s",(char*)__FILE__,__LINE__,fReplaceString1.Data(),
1c21a0e7 217 newName.Data());
fb10c4b8 218 continue;
219 }
220
ea950747 221
fb10c4b8 222 // add a new branch to the output for the background subtracted jets take the names from
223 // the input jets and replace the background flag names
224 TClonesArray *tca = new TClonesArray("AliAODJet", 0);
fb10c4b8 225 newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
ea950747 226 if(newName.CompareTo(ostr->GetString())==0){
227 Printf("%s:%d Input and output branch would have the same name, skipping: %s ",(char*)__FILE__,__LINE__,ostr->GetString().Data());
228 continue;
229 }
230
fb10c4b8 231 if(fDebug){
232 Printf("%s:%d created branch \n %s from \n %s",(char*)__FILE__,__LINE__,newName.Data(),
1c21a0e7 233 ostr->GetString().Data());
fb10c4b8 234 }
235 tca->SetName(newName.Data());
236 AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
fb10c4b8 237 }
238
239
fb10c4b8 240 if(!fHistList)fHistList = new TList();
241 fHistList->SetOwner();
750e22e8 242 PostData(1, fHistList); // post data in any case once
fb10c4b8 243
ca6d12f9 244 //
245
246 // delta pT vs. area vs. cent vs. mult
247 const Int_t nSparseBinsDelta = 4;
248 const Int_t nBinsDelta[nSparseBinsDelta] = { 241, 10, 10, 25};
249 const Double_t xminDelta[nSparseBinsDelta] = {-120.5, 0, 0, 0};
250 const Double_t xmaxDelta[nSparseBinsDelta] = { 120.5, 1.0, 100,5000};
251
f0659f11 252 for(int iJB = 0;iJB<(fJBArray?fJBArray->GetEntries():0);iJB++){
1c21a0e7 253 TObjString *ostr = (TObjString*)fJBArray->At(iJB);
254 TString oldName(ostr->GetString().Data());
255 TString newName(ostr->GetString().Data());
256 if(!newName.Contains(fReplaceString1.Data())){
257 Printf("%s:%d cannot replace string %s in %s",(char*)__FILE__,__LINE__,fReplaceString1.Data(),
258 newName.Data());
259 continue;
260 }
261 newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
ca6d12f9 262
8a320ab9 263 TH2F *hTmp = new TH2F(Form("h2PtInPtOut_%d",iJB),Form(";%s p_{T}; %s p_{T}",oldName.Data(),newName.Data()),200,0,200.,400,-200.,200.);
d9c0c425 264 fHistList->Add(hTmp);
ca6d12f9 265 THnSparseF *hFTmp = new THnSparseF(Form("hnDPtAreaCentMult_%d",iJB),Form("%s delta;#delta p_{T};Area;cent;mult",newName.Data()),nSparseBinsDelta,nBinsDelta,xminDelta,xmaxDelta);
266 fHistList->Add(hFTmp);
d9c0c425 267 }
268
fb10c4b8 269 Bool_t oldStatus = TH1::AddDirectoryStatus();
270 TH1::AddDirectory(kFALSE);
271
272 //
273 // Histogram booking, add som control histograms here
274 //
275
d73e1768 276
ca6d12f9 277 fh2CentvsRho = new TH2F("fh2CentvsRho","centrality vs background density", 100,0.,100.,600,0.,300.);
278 fh2CentvsSigma = new TH2F("fh2CentvsSigma","centrality vs backgroun sigma",100,0.,100.,500,0.,50.);
d73e1768 279 fHistList->Add(fh2CentvsRho);
280 fHistList->Add(fh2CentvsSigma);
281
ca6d12f9 282 fh2MultvsRho = new TH2F("fh2MultvsRho","mult vs background density", 100,0.,5000.,600,0.,300.);
283 fh2MultvsSigma = new TH2F("fh2MultvsSigma","mult vs background sigma",100,0.,5000.,500,0.,50.);
82ac956f 284 fHistList->Add(fh2MultvsRho);
285 fHistList->Add(fh2MultvsSigma);
286
d73e1768 287 if(fSubtraction==k4Area){
288 fh2ShiftEta = new TH2F("fh2ShiftEta","extended correction Eta",100,-0.9,0.9,100,-0.9,0.9);
289 fh2ShiftPhi = new TH2F("fh2ShiftPhi","extended correction Phi",100,0.,6.5,100,0.,6.5);
290 fh2ShiftEtaLeading = new TH2F("fh2ShiftEtaLeading","extended correction Eta",100,-0.9,0.9,100,-0.9,0.9);
291 fh2ShiftPhiLeading = new TH2F("fh2ShiftPhiLeading","extended correction Phi",100,0.,6.5,100,0.,6.5);
292 fHistList->Add(fh2ShiftEta);
293 fHistList->Add(fh2ShiftPhi);
294 fHistList->Add(fh2ShiftEtaLeading);
295 fHistList->Add(fh2ShiftPhiLeading);
296 }
297
fb10c4b8 298 // =========== Switch on Sumw2 for all histos ===========
299 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
300 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
301 if (h1){
302 h1->Sumw2();
303 continue;
304 }
305 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
306 if(hn)hn->Sumw2();
307 }
308 TH1::AddDirectory(oldStatus);
309
310 if(fBackgroundBranch.Length()==0)
311 AliError(Form("%s:%d No BackgroundBranch defined",(char*)__FILE__,__LINE__));
f0659f11 312 if((fJBArray?fJBArray->GetEntries():0)==0)
fb10c4b8 313 AliError(Form("%s:%d No Jet Branches defined defined",(char*)__FILE__,__LINE__));
314}
315
316void AliAnalysisTaskJetBackgroundSubtract::Init()
317{
318 //
319 // Initialization
320 //
321 if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::Init() \n");
322}
323
324void AliAnalysisTaskJetBackgroundSubtract::UserExec(Option_t */*option*/)
325{
f0659f11 326 //
327 // Execute for every selected event
328 //
fb10c4b8 329
1c21a0e7 330 if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserExec() \n");
fb10c4b8 331 ResetOutJets();
f0659f11 332 if(fBackgroundBranch.Length()==0||(fJBArray?fJBArray->GetEntries():0)==0){
1c21a0e7 333 if(fDebug)Printf("%s:%d No background subtraction done",(char*)__FILE__,__LINE__);
fb10c4b8 334 PostData(1,fHistList);
335 }
1c21a0e7 336 if(fJBArray->GetEntries()!=fInJetArrayList->GetEntries()){
337 if(fDebug)Printf("%s:%d different Array sizes %d %d %d",(char*)__FILE__,__LINE__,fJBArray->GetEntries(),fInJetArrayList->GetEntries(),fOutJetArrayList->GetEntries());
338 PostData(1,fHistList);
339 }
340
fb10c4b8 341
342
a18eedbb 343 AliAODJetEventBackground* evBkg = 0;
344 TClonesArray* bkgClusters = 0;
345 TClonesArray* bkgClustersRC = 0;
346 TString bkgClusterName(fBackgroundBranch.Data());
1c21a0e7 347 bkgClusterName.ReplaceAll(Form("%s_",AliAODJetEventBackground::StdBranchName()),"");
a18eedbb 348 TString bkgClusterRCName(Form("%s%s",bkgClusterName.Data(),"RandomCone"));
349
350 if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODOut){
fb10c4b8 351 evBkg = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
d9c0c425 352 bkgClusters = (TClonesArray*)(fAODOut->FindListObject(bkgClusterName.Data()));
a18eedbb 353 bkgClustersRC = (TClonesArray*)(fAODOut->FindListObject(bkgClusterRCName.Data()));
354
1c21a0e7 355 if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
a18eedbb 356 if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
1c21a0e7 357 if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
fb10c4b8 358 }
a18eedbb 359 if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODExtension){
fb10c4b8 360 evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
d9c0c425 361 bkgClusters = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(bkgClusterName.Data()));
a18eedbb 362 bkgClustersRC = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(bkgClusterRCName.Data()));
1c21a0e7 363 if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
a18eedbb 364 if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
365
1c21a0e7 366 if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
fb10c4b8 367 }
1c21a0e7 368
a18eedbb 369 if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODIn){
fb10c4b8 370 evBkg = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
1c21a0e7 371 bkgClusters = (TClonesArray*)(fAODIn->FindListObject(bkgClusterName.Data()));
a18eedbb 372 bkgClustersRC = (TClonesArray*)(fAODIn->FindListObject(bkgClusterRCName.Data()));
373
1c21a0e7 374 if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
a18eedbb 375 if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
1c21a0e7 376 if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
fb10c4b8 377 }
378
8ddcf605 379 if(!evBkg&&(fSubtraction==kArea||fSubtraction==kRhoRecalc||fSubtraction==k4Area)){
1c21a0e7 380 if(fDebug){
381 Printf("%s:%d Backgroundbranch %s not found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
382 PrintAODContents();
383 }
fb10c4b8 384 PostData(1,fHistList);
d9c0c425 385 return;
386 }
387
8a56597f 388 if(!bkgClusters&&(fSubtraction==kRhoRecalc)){
1c21a0e7 389 if(fDebug){
390 Printf("%s:%d Background cluster branch %s not found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
391 PrintAODContents();
392 }
d9c0c425 393 PostData(1,fHistList);
394 return;
fb10c4b8 395 }
396
8a56597f 397 if(!bkgClustersRC&&(fSubtraction==kRhoRC)){
a18eedbb 398 if(fDebug){
399 Printf("%s:%d Background cluster RC branch %s not found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
400 PrintAODContents();
401 }
402 PostData(1,fHistList);
403 return;
404 }
fb10c4b8 405 // LOOP over all jet branches and subtract the background
406
a18eedbb 407 Float_t rho = 0;
d73e1768 408 Float_t sigma=0.;
a18eedbb 409 Double_t meanarea = 0;
d73e1768 410 TLorentzVector backgroundv;
411 Float_t cent=0.;
82ac956f 412
8ddcf605 413 if(fAODOut)cent = fAODOut->GetHeader()->GetCentrality();
824a37c1 414 if(fAODIn) cent = fAODIn->GetHeader()->GetCentrality();
415
8ddcf605 416 if(evBkg)sigma=evBkg->GetSigma(1);
e383bb7b 417
418 if(fSubtraction==kArea) {rho = evBkg->GetBackground(1);
419 if(fExclude2Leading==kFALSE){
420 rho=evBkg->GetBackground(3);
421 sigma=evBkg->GetSigma(3);}}
d73e1768 422
d73e1768 423 if(fSubtraction==k4Area){
e383bb7b 424
8ddcf605 425 rho = evBkg->GetBackground(0);
e383bb7b 426 sigma=evBkg->GetSigma(0);
8ddcf605 427 }
e383bb7b 428
429
d73e1768 430
a18eedbb 431 if(fSubtraction==kRhoRecalc){
d73e1768 432 meanarea=evBkg->GetMeanarea(1);
a18eedbb 433 rho =RecalcRho(bkgClusters,meanarea);
8ddcf605 434 }
d73e1768 435 if(fSubtraction==kRhoRC) rho=RhoRC(bkgClustersRC);
82ac956f 436
437 Float_t mult = 0;
438 for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
439 TClonesArray* jarray = (TClonesArray*)fInJetArrayList->At(iJB);
440 if(jarray){
441 TString tmp(jarray->GetName());
442 if(tmp.Contains("cluster")){
443 mult = MultFromJetRefs(jarray);
444 if(mult>0)break;
445 }
446 }
447 }
448
d73e1768 449 fh2CentvsRho->Fill(cent,rho);
450 fh2CentvsSigma->Fill(cent,sigma);
82ac956f 451
452 fh2MultvsRho->Fill(mult,rho);
453 fh2MultvsSigma->Fill(mult,sigma);
a18eedbb 454
d73e1768 455 for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
fb10c4b8 456 TClonesArray* jarray = (TClonesArray*)fInJetArrayList->At(iJB);
457 TClonesArray* jarrayOut = (TClonesArray*)fOutJetArrayList->At(iJB);
1c21a0e7 458
459 if(!jarray||!jarrayOut){
8ddcf605 460 Printf("%s:%d Array not found %d: %p %p",(char*)__FILE__,__LINE__,iJB,jarray,jarrayOut);
461 continue;
462 }
d9c0c425 463 TH2F* h2PtInOut = (TH2F*)fHistList->FindObject(Form("h2PtInPtOut_%d",iJB));
ca6d12f9 464 THnSparseF* hnDPtAreaCentMult = (THnSparseF*)fHistList->FindObject(Form("hnDPtAreaCentMult_%d",iJB));
fb10c4b8 465 // loop over all jets
466 Int_t nOut = 0;
a18eedbb 467
ca6d12f9 468 Double_t deltaPt[4];
469 deltaPt[2] = cent;
470 deltaPt[3] = mult;
a18eedbb 471
fb10c4b8 472 for(int i = 0;i < jarray->GetEntriesFast();i++){
473 AliAODJet *jet = (AliAODJet*)jarray->At(i);
474 AliAODJet tmpNewJet(*jet);
475 Bool_t bAdd = false;
ca6d12f9 476 Float_t ptSub = 0;
477
a18eedbb 478
fb10c4b8 479 if(fSubtraction==kArea){
d9c0c425 480 Double_t background = rho * jet->EffectiveAreaCharged();
ca6d12f9 481 ptSub = jet->Pt() - background;
fb10c4b8 482 if(fDebug>2){
483 Printf("%s:%d Jet %d %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub);
484 }
c0711eb9 485 if(ptSub<=0){
fb10c4b8 486 // optionally rescale it and keep??
bea24d2d 487 if(fKeepJets){
488 bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
489 }
490 else{
491 bAdd = false;
492 }
fb10c4b8 493 }
494 else{
495 bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
496 }
d9c0c425 497 // add background estimates to the new jet object
498 // allows to recover old p_T and rho...
499 tmpNewJet.SetBgEnergy(background,0);
bea24d2d 500 tmpNewJet.SetPtSubtracted(ptSub,0);
fb10c4b8 501 }// kAREA
a18eedbb 502 else if(fSubtraction==kRhoRecalc){
503 Double_t background = rho * jet->EffectiveAreaCharged();
ca6d12f9 504 ptSub = jet->Pt() - background;
a18eedbb 505 if(fDebug>2){
506 Printf("%s:%d Jet %d %3.3f %3.3f %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub,background,rho);}
c0711eb9 507 if(ptSub<=0){
bea24d2d 508 // optionally rescale it and keep
509 if(fKeepJets){
510 bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
511 }
512 else{
513 bAdd = false;
514 }
a18eedbb 515 }
516 else{
517 bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
a18eedbb 518 }
519 // add background estimates to the new jet object
520 // allows to recover old p_T and rho...
521 tmpNewJet.SetBgEnergy(background,0);
bea24d2d 522 tmpNewJet.SetPtSubtracted(ptSub,0);
a18eedbb 523 }//kRhoRecalc
524 else if(fSubtraction==kRhoRC){
525 Double_t background = rho * jet->EffectiveAreaCharged();
ca6d12f9 526 ptSub = jet->Pt() - background;
a18eedbb 527 if(fDebug>2){ Printf("%s:%d Jet %d %3.3f %3.3f %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub,background,rho);}
c0711eb9 528 if(ptSub<=0){
bea24d2d 529 if(fKeepJets){
530 bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
531 }
532 else{
533 bAdd = false;
534 }
a18eedbb 535 }
536 else{
537 bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
a18eedbb 538 }
539 // add background estimates to the new jet object
540 // allows to recover old p_T and rho...
541 tmpNewJet.SetBgEnergy(background,0);
bea24d2d 542 tmpNewJet.SetPtSubtracted(ptSub,0);
a18eedbb 543 }//kRhoRC
fb10c4b8 544
c0ec1c44 545 else if(fSubtraction==k4Area&&jet->VectorAreaCharged()){
c0ec1c44 546 backgroundv.SetPxPyPzE(rho*(jet->VectorAreaCharged())->Px(),rho*(jet->VectorAreaCharged())->Py(),rho*(jet->VectorAreaCharged())->Pz(),rho*(jet->VectorAreaCharged())->E());
ca6d12f9 547 ptSub = jet->Pt()-backgroundv.Pt();
c0711eb9 548 if((backgroundv.E()>=jet->E())||(backgroundv.Pt()>=jet->Pt())){
bea24d2d 549 if(fKeepJets){
550 bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
551 }
552 else{
553 bAdd = false;
554 }
c0ec1c44 555 }
556 else{
557 bAdd = RescaleJet4vector(&tmpNewJet,backgroundv);
558 }
559 // add background estimates to the new jet object
560 // allows to recover old p_T and rho...
bea24d2d 561 tmpNewJet.SetBgEnergy(backgroundv.Pt(),0);
562 tmpNewJet.SetPtSubtracted(ptSub,0);
c0ec1c44 563
d73e1768 564 }//kArea4vector
565
fb10c4b8 566 if(bAdd){
d9c0c425 567 AliAODJet *newJet = new ((*jarrayOut)[nOut++]) AliAODJet(tmpNewJet);
568 // what about track references, clear for now...
d73e1768 569 if(fSubtraction==k4Area){
d73e1768 570 fh2ShiftEta->Fill(jet->Eta(),newJet->Eta());
571 fh2ShiftPhi->Fill(jet->Phi(),newJet->Phi());
572 if(i==0){fh2ShiftEtaLeading->Fill(jet->Eta(),newJet->Eta());
573 fh2ShiftPhiLeading->Fill(jet->Phi(),newJet->Phi());}}
574
d508795b 575 // set the references
d9c0c425 576 newJet->GetRefTracks()->Clear();
d508795b 577 TRefArray *refs = jet->GetRefTracks();
578 for(Int_t ir=0;ir<refs->GetEntriesFast();ir++){
579 AliVParticle *vp = dynamic_cast<AliVParticle*>(refs->At(ir));
580 if(vp)newJet->AddTrack(vp);
581 }
fb10c4b8 582 }
ca6d12f9 583 if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
584 if(hnDPtAreaCentMult){
585 deltaPt[0] = ptSub;
586 deltaPt[1] = jet->EffectiveAreaCharged();
587 hnDPtAreaCentMult->Fill(deltaPt);
588 }
a18eedbb 589 }
cb9ef56f 590 if(jarrayOut)jarrayOut->Sort();
591 }
592
593 PostData(1, fHistList);
fb10c4b8 594}
595
596void AliAnalysisTaskJetBackgroundSubtract::Terminate(Option_t */*option*/)
597{
598 // Terminate analysis
599 //
600 if (fDebug > 1) printf("AnalysisJetBackgroundSubtract: Terminate() \n");
601}
602
603Bool_t AliAnalysisTaskJetBackgroundSubtract::RescaleJetMomentum(AliAODJet *jet,Float_t pT){
604 // keep the direction and the jet mass
605 if(pT<=0)return kFALSE;
606 Double_t pTold = jet->Pt();
607 Double_t scale = pT/pTold;
608 Double_t mass = jet->M();
609 Double_t pNew = jet->P() * scale;
610 jet->SetPxPyPzE(scale*jet->Px(),scale*jet->Py(),scale*jet->Pz(),TMath::Sqrt(mass*mass+pNew*pNew));
d73e1768 611
612
613
fb10c4b8 614 return kTRUE;
615}
616
d73e1768 617Bool_t AliAnalysisTaskJetBackgroundSubtract::RescaleJet4vector(AliAODJet *jet,TLorentzVector backgroundv){
618
619 if(backgroundv.Pt()<0.) return kFALSE;
620 jet->SetPxPyPzE(jet->Px()-backgroundv.Px(),jet->Py()-backgroundv.Py(),jet->Pz()-backgroundv.Pz(),jet->E()-backgroundv.E());
621
622 return kTRUE;
623}
624
625
626
627
628
629
630
631
a18eedbb 632Double_t AliAnalysisTaskJetBackgroundSubtract::RecalcRho(TClonesArray* bkgClusters,Double_t meanarea){
633
f0659f11 634 //
635 // recalc rhoo
636 //
637
8ddcf605 638 Double_t ptarea=0.;
639 Int_t count=0;
640 Double_t rho=0.;
f0659f11 641 const Double_t rLimit2=0.8*0.8; //2*jet radius.
8ddcf605 642 TClonesArray* jarray=0;
643
644 for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
645 TObjString *ostr = (TObjString*)fInJetArrayList->At(iJB);
646 TString jetref=ostr->GetString().Data();
647 if(jetref.Contains("ANTIKT04")){
648 jarray = (TClonesArray*)fInJetArrayList->At(iJB);
649 }
650 }
651 if(!jarray)return rho;
652 if(jarray->GetEntries()>=2){
653 AliAODJet *first = (AliAODJet*)(jarray->At(0));
654 AliAODJet *second= (AliAODJet*)(jarray->At(1));
655 for(Int_t k=0;k<bkgClusters->GetEntriesFast();k++){
656 AliAODJet *clus = (AliAODJet*)(bkgClusters->At(k));
657 if(TMath::Abs(clus->Eta())>0.5) continue;
658 if((clus->EffectiveAreaCharged())<0.1*meanarea) continue;
659 Double_t distance1=(first->Eta()-clus->Eta())*(first->Eta()-clus->Eta())+
660 (first->Phi()-clus->Phi())*(first->Phi()-clus->Phi());
661 Double_t distance2= (second->Eta()-clus->Eta())*(second->Eta()-clus->Eta())+
662 (second->Phi()-clus->Phi())*(second->Phi()-clus->Phi());
f0659f11 663 if((distance1<rLimit2)||(distance2<rLimit2)) continue;
8ddcf605 664 ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged();
665 count=count+1;}
666 if(count!=0) rho=ptarea/count;
667 }
668 return rho;
a18eedbb 669}
670
f0659f11 671Double_t AliAnalysisTaskJetBackgroundSubtract::RhoRC(TClonesArray* bkgClustersRC){
672
673 //
674 // calc rho from random cones
675 //
a18eedbb 676
677 Double_t ptarea=0.;
678 Int_t count=0;
679 Double_t rho=0.;
f0659f11 680 const Double_t rLimit2=0.8*0.8; //2*jet radius.
a18eedbb 681 TClonesArray* jarray=0;
682 for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
8ddcf605 683 TObjString *ostr = (TObjString*)fInJetArrayList->At(iJB);
684 TString jetref=ostr->GetString().Data();
685 if(jetref.Contains("ANTIKT04")){
686 jarray = (TClonesArray*)fInJetArrayList->At(iJB);
687 }
688 }
689 if(!jarray)return rho;
690
a18eedbb 691 if(jarray->GetEntries()>=2){
692 AliAODJet *first = (AliAODJet*)(jarray->At(0));
693 AliAODJet *second=(AliAODJet*)(jarray->At(1));
694 for(Int_t k=0;k<bkgClustersRC->GetEntriesFast();k++){
695 AliAODJet *clus = (AliAODJet*)(bkgClustersRC->At(k));
696 if(TMath::Abs(clus->Eta())>0.5) continue;
697 Double_t distance1=(first->Eta()-clus->Eta())*(first->Eta()-clus->Eta())+
698 (first->Phi()-clus->Phi())*(first->Phi()-clus->Phi());
699 Double_t distance2= (second->Eta()-clus->Eta())*(second->Eta()-clus->Eta())+
700 (second->Phi()-clus->Phi())*(second->Phi()-clus->Phi());
f0659f11 701 if((distance1<rLimit2)||(distance2<rLimit2)) continue;
a18eedbb 702 ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged();
703 count=count+1;}
704 if(count!=0) rho=ptarea/count; }
705 return rho;
706}
707
708
709
710
711
712
713
714
715
fb10c4b8 716void AliAnalysisTaskJetBackgroundSubtract::ResetOutJets(){
f0659f11 717 //
718 // Reset the output jets
719 //
720
fb10c4b8 721 if(!fOutJetArrayList)return;
722 for(int iJB = 0;iJB<fOutJetArrayList->GetEntries();iJB++){
723 TClonesArray* jarray = (TClonesArray*)fOutJetArrayList->At(iJB);
724 if(jarray)jarray->Delete();
725 }
726}
1c21a0e7 727
728
729void AliAnalysisTaskJetBackgroundSubtract::PrintAODContents(){
f0659f11 730
731 //
732 // guess from the name what this function does
733 //
734
1c21a0e7 735 if(fAODIn){
736 Printf("%s:%d >>>>>> Input",(char*)__FILE__,__LINE__);
737 fAODIn->Print();
738 }
739 if(fAODExtension){
740 Printf("%s:%d >>>>>> Extenstion",(char*)__FILE__,__LINE__);
741 fAODExtension->GetAOD()->Print();
742 }
743 if(fAODOut){
744 Printf("%s:%d >>>>>> Output",(char*)__FILE__,__LINE__);
745 fAODOut->Print();
746 }
747}
82ac956f 748
749Int_t AliAnalysisTaskJetBackgroundSubtract::MultFromJetRefs(TClonesArray* jets){
f0659f11 750 //
751 // calculate multiplicty based on jet references
752 //
753
82ac956f 754 if(!jets)return 0;
755
756 Int_t refMult = 0;
757 for(int ij = 0;ij < jets->GetEntries();++ij){
758 AliAODJet* jet = (AliAODJet*)jets->At(ij);
759 if(!jet)continue;
760 TRefArray *refs = jet->GetRefTracks();
761 if(!refs)continue;
762 refMult += refs->GetEntries();
763 }
764 return refMult;
765
766}
f0659f11 767
768void AliAnalysisTaskJetBackgroundSubtract::AddJetBranch(const char* c){
769 if(!fJBArray)fJBArray = new TObjArray();
770 fJBArray->Add(new TObjString(c));
771}