]> 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
adapted histo ranges
[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();
380 if(evBkg)sigma=evBkg->GetSigma(1);
d73e1768 381
382 if(fSubtraction==kArea) rho = evBkg->GetBackground(1);
383 if(fSubtraction==k4Area){
8ddcf605 384 rho = evBkg->GetBackground(0);
385 sigma=evBkg->GetSigma(0);
386 }
d73e1768 387
a18eedbb 388 if(fSubtraction==kRhoRecalc){
d73e1768 389 meanarea=evBkg->GetMeanarea(1);
a18eedbb 390 rho =RecalcRho(bkgClusters,meanarea);
8ddcf605 391 }
d73e1768 392 if(fSubtraction==kRhoRC) rho=RhoRC(bkgClustersRC);
82ac956f 393
394 Float_t mult = 0;
395 for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
396 TClonesArray* jarray = (TClonesArray*)fInJetArrayList->At(iJB);
397 if(jarray){
398 TString tmp(jarray->GetName());
399 if(tmp.Contains("cluster")){
400 mult = MultFromJetRefs(jarray);
401 if(mult>0)break;
402 }
403 }
404 }
405
d73e1768 406 fh2CentvsRho->Fill(cent,rho);
407 fh2CentvsSigma->Fill(cent,sigma);
82ac956f 408
409 fh2MultvsRho->Fill(mult,rho);
410 fh2MultvsSigma->Fill(mult,sigma);
a18eedbb 411
d73e1768 412 for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
fb10c4b8 413 TClonesArray* jarray = (TClonesArray*)fInJetArrayList->At(iJB);
414 TClonesArray* jarrayOut = (TClonesArray*)fOutJetArrayList->At(iJB);
1c21a0e7 415
416 if(!jarray||!jarrayOut){
8ddcf605 417 Printf("%s:%d Array not found %d: %p %p",(char*)__FILE__,__LINE__,iJB,jarray,jarrayOut);
418 continue;
419 }
d9c0c425 420 TH2F* h2PtInOut = (TH2F*)fHistList->FindObject(Form("h2PtInPtOut_%d",iJB));
fb10c4b8 421 // loop over all jets
422 Int_t nOut = 0;
a18eedbb 423
424
fb10c4b8 425 for(int i = 0;i < jarray->GetEntriesFast();i++){
426 AliAODJet *jet = (AliAODJet*)jarray->At(i);
427 AliAODJet tmpNewJet(*jet);
428 Bool_t bAdd = false;
d73e1768 429
a18eedbb 430
fb10c4b8 431 if(fSubtraction==kArea){
d9c0c425 432 Double_t background = rho * jet->EffectiveAreaCharged();
433 Float_t ptSub = jet->Pt() - background;
fb10c4b8 434 if(fDebug>2){
435 Printf("%s:%d Jet %d %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub);
436 }
437 if(ptSub<0){
438 // optionally rescale it and keep??
8a320ab9 439 bAdd = false; // RescaleJetMomentum(&tmpNewJet,0.1);
440 if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
fb10c4b8 441 }
442 else{
443 bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
d9c0c425 444 if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
fb10c4b8 445 }
d9c0c425 446 // add background estimates to the new jet object
447 // allows to recover old p_T and rho...
448 tmpNewJet.SetBgEnergy(background,0);
fb10c4b8 449 }// kAREA
a18eedbb 450 else if(fSubtraction==kRhoRecalc){
451 Double_t background = rho * jet->EffectiveAreaCharged();
452 Float_t ptSub = jet->Pt() - background;
453 if(fDebug>2){
454 Printf("%s:%d Jet %d %3.3f %3.3f %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub,background,rho);}
455 if(ptSub<0){
456 // optionally rescale it and keep??
8a320ab9 457 bAdd = false;// RescaleJetMomentum(&tmpNewJet,0.1);
458 if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
a18eedbb 459 }
460 else{
461 bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
462 if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
463 }
464 // add background estimates to the new jet object
465 // allows to recover old p_T and rho...
466 tmpNewJet.SetBgEnergy(background,0);
467
468 }//kRhoRecalc
469 else if(fSubtraction==kRhoRC){
470 Double_t background = rho * jet->EffectiveAreaCharged();
471 Float_t ptSub = jet->Pt() - background;
472 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);}
473 if(ptSub<0){
474 // optionally rescale it and keep??
8a320ab9 475 bAdd = false; // RescaleJetMomentum(&tmpNewJet,0.1);
476 if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
a18eedbb 477 }
478 else{
479 bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
480 if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
481 }
482 // add background estimates to the new jet object
483 // allows to recover old p_T and rho...
484 tmpNewJet.SetBgEnergy(background,0);
485
486 }//kRhoRC
fb10c4b8 487
c0ec1c44 488 else if(fSubtraction==k4Area&&jet->VectorAreaCharged()){
d73e1768 489
c0ec1c44 490
491 backgroundv.SetPxPyPzE(rho*(jet->VectorAreaCharged())->Px(),rho*(jet->VectorAreaCharged())->Py(),rho*(jet->VectorAreaCharged())->Pz(),rho*(jet->VectorAreaCharged())->E());
492 if((backgroundv.E()>jet->E())&&(backgroundv.Pt()>jet->Pt())){
493
494 // optionally rescale it and keep??
8a320ab9 495 bAdd = false; // RescaleJetMomentum(&tmpNewJet,0.1);
496 if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),jet->Pt()-backgroundv.Pt());
c0ec1c44 497 }
498 else{
499 bAdd = RescaleJet4vector(&tmpNewJet,backgroundv);
500 }
501 // add background estimates to the new jet object
502 // allows to recover old p_T and rho...
503 tmpNewJet.SetBgEnergy(backgroundv.P(),0);
504
d73e1768 505 }//kArea4vector
506
507
508
509
510
511
fb10c4b8 512 if(bAdd){
d9c0c425 513 AliAODJet *newJet = new ((*jarrayOut)[nOut++]) AliAODJet(tmpNewJet);
514 // what about track references, clear for now...
d73e1768 515 if(fSubtraction==k4Area){
516 if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),jet->Pt()-newJet->Pt());
517 fh2ShiftEta->Fill(jet->Eta(),newJet->Eta());
518 fh2ShiftPhi->Fill(jet->Phi(),newJet->Phi());
519 if(i==0){fh2ShiftEtaLeading->Fill(jet->Eta(),newJet->Eta());
520 fh2ShiftPhiLeading->Fill(jet->Phi(),newJet->Phi());}}
521
d508795b 522 // set the references
d9c0c425 523 newJet->GetRefTracks()->Clear();
d508795b 524 TRefArray *refs = jet->GetRefTracks();
525 for(Int_t ir=0;ir<refs->GetEntriesFast();ir++){
526 AliVParticle *vp = dynamic_cast<AliVParticle*>(refs->At(ir));
527 if(vp)newJet->AddTrack(vp);
528 }
fb10c4b8 529 }
a18eedbb 530 }
cb9ef56f 531 if(jarrayOut)jarrayOut->Sort();
532 }
533
534 PostData(1, fHistList);
fb10c4b8 535}
536
537void AliAnalysisTaskJetBackgroundSubtract::Terminate(Option_t */*option*/)
538{
539 // Terminate analysis
540 //
541 if (fDebug > 1) printf("AnalysisJetBackgroundSubtract: Terminate() \n");
542}
543
544Bool_t AliAnalysisTaskJetBackgroundSubtract::RescaleJetMomentum(AliAODJet *jet,Float_t pT){
545 // keep the direction and the jet mass
546 if(pT<=0)return kFALSE;
547 Double_t pTold = jet->Pt();
548 Double_t scale = pT/pTold;
549 Double_t mass = jet->M();
550 Double_t pNew = jet->P() * scale;
551 jet->SetPxPyPzE(scale*jet->Px(),scale*jet->Py(),scale*jet->Pz(),TMath::Sqrt(mass*mass+pNew*pNew));
d73e1768 552
553
554
fb10c4b8 555 return kTRUE;
556}
557
d73e1768 558Bool_t AliAnalysisTaskJetBackgroundSubtract::RescaleJet4vector(AliAODJet *jet,TLorentzVector backgroundv){
559
560 if(backgroundv.Pt()<0.) return kFALSE;
561 jet->SetPxPyPzE(jet->Px()-backgroundv.Px(),jet->Py()-backgroundv.Py(),jet->Pz()-backgroundv.Pz(),jet->E()-backgroundv.E());
562
563 return kTRUE;
564}
565
566
567
568
569
570
571
572
a18eedbb 573Double_t AliAnalysisTaskJetBackgroundSubtract::RecalcRho(TClonesArray* bkgClusters,Double_t meanarea){
574
8ddcf605 575 Double_t ptarea=0.;
576 Int_t count=0;
577 Double_t rho=0.;
578 const Double_t Rlimit2=0.8*0.8; //2*jet radius.
579 TClonesArray* jarray=0;
580
581 for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
582 TObjString *ostr = (TObjString*)fInJetArrayList->At(iJB);
583 TString jetref=ostr->GetString().Data();
584 if(jetref.Contains("ANTIKT04")){
585 jarray = (TClonesArray*)fInJetArrayList->At(iJB);
586 }
587 }
588 if(!jarray)return rho;
589 if(jarray->GetEntries()>=2){
590 AliAODJet *first = (AliAODJet*)(jarray->At(0));
591 AliAODJet *second= (AliAODJet*)(jarray->At(1));
592 for(Int_t k=0;k<bkgClusters->GetEntriesFast();k++){
593 AliAODJet *clus = (AliAODJet*)(bkgClusters->At(k));
594 if(TMath::Abs(clus->Eta())>0.5) continue;
595 if((clus->EffectiveAreaCharged())<0.1*meanarea) continue;
596 Double_t distance1=(first->Eta()-clus->Eta())*(first->Eta()-clus->Eta())+
597 (first->Phi()-clus->Phi())*(first->Phi()-clus->Phi());
598 Double_t distance2= (second->Eta()-clus->Eta())*(second->Eta()-clus->Eta())+
599 (second->Phi()-clus->Phi())*(second->Phi()-clus->Phi());
600 if((distance1<Rlimit2)||(distance2<Rlimit2)) continue;
601 ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged();
602 count=count+1;}
603 if(count!=0) rho=ptarea/count;
604 }
605 return rho;
a18eedbb 606}
607
608 Double_t AliAnalysisTaskJetBackgroundSubtract::RhoRC(TClonesArray* bkgClustersRC){
609
610 Double_t ptarea=0.;
611 Int_t count=0;
612 Double_t rho=0.;
613 const Double_t Rlimit2=0.8*0.8; //2*jet radius.
614 TClonesArray* jarray=0;
615 for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
8ddcf605 616 TObjString *ostr = (TObjString*)fInJetArrayList->At(iJB);
617 TString jetref=ostr->GetString().Data();
618 if(jetref.Contains("ANTIKT04")){
619 jarray = (TClonesArray*)fInJetArrayList->At(iJB);
620 }
621 }
622 if(!jarray)return rho;
623
a18eedbb 624 if(jarray->GetEntries()>=2){
625 AliAODJet *first = (AliAODJet*)(jarray->At(0));
626 AliAODJet *second=(AliAODJet*)(jarray->At(1));
627 for(Int_t k=0;k<bkgClustersRC->GetEntriesFast();k++){
628 AliAODJet *clus = (AliAODJet*)(bkgClustersRC->At(k));
629 if(TMath::Abs(clus->Eta())>0.5) continue;
630 Double_t distance1=(first->Eta()-clus->Eta())*(first->Eta()-clus->Eta())+
631 (first->Phi()-clus->Phi())*(first->Phi()-clus->Phi());
632 Double_t distance2= (second->Eta()-clus->Eta())*(second->Eta()-clus->Eta())+
633 (second->Phi()-clus->Phi())*(second->Phi()-clus->Phi());
634 if((distance1<Rlimit2)||(distance2<Rlimit2)) continue;
635 ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged();
636 count=count+1;}
637 if(count!=0) rho=ptarea/count; }
638 return rho;
639}
640
641
642
643
644
645
646
647
648
fb10c4b8 649void AliAnalysisTaskJetBackgroundSubtract::ResetOutJets(){
650 if(!fOutJetArrayList)return;
651 for(int iJB = 0;iJB<fOutJetArrayList->GetEntries();iJB++){
652 TClonesArray* jarray = (TClonesArray*)fOutJetArrayList->At(iJB);
653 if(jarray)jarray->Delete();
654 }
655}
1c21a0e7 656
657
658void AliAnalysisTaskJetBackgroundSubtract::PrintAODContents(){
659 if(fAODIn){
660 Printf("%s:%d >>>>>> Input",(char*)__FILE__,__LINE__);
661 fAODIn->Print();
662 }
663 if(fAODExtension){
664 Printf("%s:%d >>>>>> Extenstion",(char*)__FILE__,__LINE__);
665 fAODExtension->GetAOD()->Print();
666 }
667 if(fAODOut){
668 Printf("%s:%d >>>>>> Output",(char*)__FILE__,__LINE__);
669 fAODOut->Print();
670 }
671}
82ac956f 672
673Int_t AliAnalysisTaskJetBackgroundSubtract::MultFromJetRefs(TClonesArray* jets){
674 if(!jets)return 0;
675
676 Int_t refMult = 0;
677 for(int ij = 0;ij < jets->GetEntries();++ij){
678 AliAODJet* jet = (AliAODJet*)jets->At(ij);
679 if(!jet)continue;
680 TRefArray *refs = jet->GetRefTracks();
681 if(!refs)continue;
682 refMult += refs->GetEntries();
683 }
684 return refMult;
685
686}