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