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