]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - JETAN/AliAnalysisTaskJetBackgroundSubtract.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Fix: increase classdef
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskJetBackgroundSubtract.cxx
... / ...
CommitLineData
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>
25#include <TH2F.h>
26#include <THnSparse.h>
27#include <TSystem.h>
28#include <TInterpreter.h>
29#include <TList.h>
30#include <TLorentzVector.h>
31#include <TRefArray.h>
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"),
62 fSubtraction(k4Area),
63 fKeepJets(kFALSE),
64 fInJetArrayList(0x0),
65 fOutJetArrayList(0x0),
66 fh2CentvsRho(0x0),
67 fh2CentvsSigma(0x0),
68 fh2MultvsRho(0x0),
69 fh2MultvsSigma(0x0),
70 fh2ShiftEta(0x0),
71 fh2ShiftPhi(0x0),
72 fh2ShiftEtaLeading(0x0),
73 fh2ShiftPhiLeading(0x0),
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"),
90 fSubtraction(k4Area),
91 fKeepJets(kFALSE),
92 fInJetArrayList(0x0),
93 fOutJetArrayList(0x0),
94 fh2CentvsRho(0x0),
95 fh2CentvsSigma(0x0),
96 fh2MultvsRho(0x0),
97 fh2MultvsSigma(0x0),
98 fh2ShiftEta(0x0),
99 fh2ShiftPhi(0x0),
100 fh2ShiftEtaLeading(0x0),
101 fh2ShiftPhiLeading(0x0),
102 fHistList(0x0)
103{
104 DefineOutput(1, TList::Class());
105}
106
107
108
109Bool_t AliAnalysisTaskJetBackgroundSubtract::Notify()
110{
111 //
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);
123
124
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;
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
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 }
170 return kTRUE;
171}
172
173void AliAnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects()
174{
175
176 //
177 // Create the output container
178 //
179 // Connect the AOD
180
181 if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects() \n");
182
183
184
185 if(fNonStdFile.Length()!=0){
186
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());
190
191 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
192
193 if(!fAODExtension){
194 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
195 }
196 }
197 fAODOut = AODEvent();
198
199 // usually we do not have the input already here
200
201 if(!fInJetArrayList)fInJetArrayList =new TList();
202 if(!fOutJetArrayList)fOutJetArrayList =new TList();
203
204 for(int iJB = 0;iJB<fJBArray->GetEntries();iJB++){
205 TObjString *ostr = (TObjString*)fJBArray->At(iJB);
206 TString newName(ostr->GetString().Data());
207 if(!newName.Contains(fReplaceString1.Data())){
208 Printf("%s:%d cannot replace string %s in %s",(char*)__FILE__,__LINE__,fReplaceString1.Data(),
209 newName.Data());
210 continue;
211 }
212
213
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);
217 newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
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
223 if(fDebug){
224 Printf("%s:%d created branch \n %s from \n %s",(char*)__FILE__,__LINE__,newName.Data(),
225 ostr->GetString().Data());
226 }
227 tca->SetName(newName.Data());
228 AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
229 }
230
231
232 if(!fHistList)fHistList = new TList();
233 fHistList->SetOwner();
234 PostData(1, fHistList); // post data in any case once
235
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
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));
254
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.);
256 fHistList->Add(hTmp);
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);
259 }
260
261 Bool_t oldStatus = TH1::AddDirectoryStatus();
262 TH1::AddDirectory(kFALSE);
263
264 //
265 // Histogram booking, add som control histograms here
266 //
267
268
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.);
271 fHistList->Add(fh2CentvsRho);
272 fHistList->Add(fh2CentvsSigma);
273
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.);
276 fHistList->Add(fh2MultvsRho);
277 fHistList->Add(fh2MultvsSigma);
278
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
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
319 if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserExec() \n");
320 ResetOutJets();
321 if(fBackgroundBranch.Length()==0||fJBArray->GetEntries()==0){
322 if(fDebug)Printf("%s:%d No background subtraction done",(char*)__FILE__,__LINE__);
323 PostData(1,fHistList);
324 }
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
330
331
332 AliAODJetEventBackground* evBkg = 0;
333 TClonesArray* bkgClusters = 0;
334 TClonesArray* bkgClustersRC = 0;
335 TString bkgClusterName(fBackgroundBranch.Data());
336 bkgClusterName.ReplaceAll(Form("%s_",AliAODJetEventBackground::StdBranchName()),"");
337 TString bkgClusterRCName(Form("%s%s",bkgClusterName.Data(),"RandomCone"));
338
339 if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODOut){
340 evBkg = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
341 bkgClusters = (TClonesArray*)(fAODOut->FindListObject(bkgClusterName.Data()));
342 bkgClustersRC = (TClonesArray*)(fAODOut->FindListObject(bkgClusterRCName.Data()));
343
344 if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
345 if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
346 if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
347 }
348 if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODExtension){
349 evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
350 bkgClusters = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(bkgClusterName.Data()));
351 bkgClustersRC = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(bkgClusterRCName.Data()));
352 if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
353 if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
354
355 if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
356 }
357
358 if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODIn){
359 evBkg = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
360 bkgClusters = (TClonesArray*)(fAODIn->FindListObject(bkgClusterName.Data()));
361 bkgClustersRC = (TClonesArray*)(fAODIn->FindListObject(bkgClusterRCName.Data()));
362
363 if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
364 if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
365 if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
366 }
367
368 if(!evBkg&&(fSubtraction==kArea||fSubtraction==kRhoRecalc||fSubtraction==k4Area)){
369 if(fDebug){
370 Printf("%s:%d Backgroundbranch %s not found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
371 PrintAODContents();
372 }
373 PostData(1,fHistList);
374 return;
375 }
376
377 if(!bkgClusters&&(fSubtraction==kRhoRecalc)){
378 if(fDebug){
379 Printf("%s:%d Background cluster branch %s not found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
380 PrintAODContents();
381 }
382 PostData(1,fHistList);
383 return;
384 }
385
386 if(!bkgClustersRC&&(fSubtraction==kRhoRC)){
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 }
394 // LOOP over all jet branches and subtract the background
395
396 Float_t rho = 0;
397 Float_t sigma=0.;
398 Double_t meanarea = 0;
399 TLorentzVector backgroundv;
400 Float_t cent=0.;
401
402 if(fAODOut)cent = fAODOut->GetHeader()->GetCentrality();
403 if(fAODIn) cent = fAODIn->GetHeader()->GetCentrality();
404
405 if(evBkg)sigma=evBkg->GetSigma(1);
406
407 if(fSubtraction==kArea) rho = evBkg->GetBackground(1);
408 if(fSubtraction==k4Area){
409 rho = evBkg->GetBackground(0);
410 sigma=evBkg->GetSigma(0);
411 }
412
413 if(fSubtraction==kRhoRecalc){
414 meanarea=evBkg->GetMeanarea(1);
415 rho =RecalcRho(bkgClusters,meanarea);
416 }
417 if(fSubtraction==kRhoRC) rho=RhoRC(bkgClustersRC);
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
431 fh2CentvsRho->Fill(cent,rho);
432 fh2CentvsSigma->Fill(cent,sigma);
433
434 fh2MultvsRho->Fill(mult,rho);
435 fh2MultvsSigma->Fill(mult,sigma);
436
437 for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
438 TClonesArray* jarray = (TClonesArray*)fInJetArrayList->At(iJB);
439 TClonesArray* jarrayOut = (TClonesArray*)fOutJetArrayList->At(iJB);
440
441 if(!jarray||!jarrayOut){
442 Printf("%s:%d Array not found %d: %p %p",(char*)__FILE__,__LINE__,iJB,jarray,jarrayOut);
443 continue;
444 }
445 TH2F* h2PtInOut = (TH2F*)fHistList->FindObject(Form("h2PtInPtOut_%d",iJB));
446 THnSparseF* hnDPtAreaCentMult = (THnSparseF*)fHistList->FindObject(Form("hnDPtAreaCentMult_%d",iJB));
447 // loop over all jets
448 Int_t nOut = 0;
449
450 Double_t deltaPt[4];
451 deltaPt[2] = cent;
452 deltaPt[3] = mult;
453
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;
458 Float_t ptSub = 0;
459
460
461 if(fSubtraction==kArea){
462 Double_t background = rho * jet->EffectiveAreaCharged();
463 ptSub = jet->Pt() - background;
464 if(fDebug>2){
465 Printf("%s:%d Jet %d %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub);
466 }
467 if(ptSub<0){
468 // optionally rescale it and keep??
469 if(fKeepJets){
470 bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
471 }
472 else{
473 bAdd = false;
474 }
475 }
476 else{
477 bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
478 }
479 // add background estimates to the new jet object
480 // allows to recover old p_T and rho...
481 tmpNewJet.SetBgEnergy(background,0);
482 tmpNewJet.SetPtSubtracted(ptSub,0);
483 }// kAREA
484 else if(fSubtraction==kRhoRecalc){
485 Double_t background = rho * jet->EffectiveAreaCharged();
486 ptSub = jet->Pt() - background;
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);}
489 if(ptSub<0){
490 // optionally rescale it and keep
491 if(fKeepJets){
492 bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
493 }
494 else{
495 bAdd = false;
496 }
497 }
498 else{
499 bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
500 }
501 // add background estimates to the new jet object
502 // allows to recover old p_T and rho...
503 tmpNewJet.SetBgEnergy(background,0);
504 tmpNewJet.SetPtSubtracted(ptSub,0);
505 }//kRhoRecalc
506 else if(fSubtraction==kRhoRC){
507 Double_t background = rho * jet->EffectiveAreaCharged();
508 ptSub = jet->Pt() - background;
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);}
510 if(ptSub<0){
511 if(fKeepJets){
512 bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
513 }
514 else{
515 bAdd = false;
516 }
517 }
518 else{
519 bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
520 }
521 // add background estimates to the new jet object
522 // allows to recover old p_T and rho...
523 tmpNewJet.SetBgEnergy(background,0);
524 tmpNewJet.SetPtSubtracted(ptSub,0);
525 }//kRhoRC
526
527 else if(fSubtraction==k4Area&&jet->VectorAreaCharged()){
528 backgroundv.SetPxPyPzE(rho*(jet->VectorAreaCharged())->Px(),rho*(jet->VectorAreaCharged())->Py(),rho*(jet->VectorAreaCharged())->Pz(),rho*(jet->VectorAreaCharged())->E());
529 ptSub = jet->Pt()-backgroundv.Pt();
530 if((backgroundv.E()>jet->E())&&(backgroundv.Pt()>jet->Pt())){
531 if(fKeepJets){
532 bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
533 }
534 else{
535 bAdd = false;
536 }
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...
543 tmpNewJet.SetBgEnergy(backgroundv.Pt(),0);
544 tmpNewJet.SetPtSubtracted(ptSub,0);
545
546 }//kArea4vector
547
548 if(bAdd){
549 AliAODJet *newJet = new ((*jarrayOut)[nOut++]) AliAODJet(tmpNewJet);
550 // what about track references, clear for now...
551 if(fSubtraction==k4Area){
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
557 // set the references
558 newJet->GetRefTracks()->Clear();
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 }
564 }
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 }
571 }
572 if(jarrayOut)jarrayOut->Sort();
573 }
574
575 PostData(1, fHistList);
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));
593
594
595
596 return kTRUE;
597}
598
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
614Double_t AliAnalysisTaskJetBackgroundSubtract::RecalcRho(TClonesArray* bkgClusters,Double_t meanarea){
615
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;
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++){
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
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
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}
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}
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}