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