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