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