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