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