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