]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/DEV/AliAnalysisTaskJetBackgroundSubtract.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
file with geometry definitions, added temporarily until the one in EVE/alice-data...
[u/mrichter/AliRoot.git] / JETAN / DEV / AliAnalysisTaskJetBackgroundSubtract.cxx
CommitLineData
d89b8229 1// **************************************
2// Task used for the correction of determiantion of reconstructed jet spectra
3// Compares input (gen) and output (rec) jets
4// *******************************************
5
6
7/**************************************************************************
8 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9 * *
10 * Author: The ALICE Off-line Project. *
11 * Contributors are mentioned in the code where appropriate. *
12 * *
13 * Permission to use, copy, modify and distribute this software and its *
14 * documentation strictly for non-commercial purposes is hereby granted *
15 * without fee, provided that the above copyright notice appears in all *
16 * copies and that both the copyright notice and this permission notice *
17 * appear in the supporting documentation. The authors make no claims *
18 * about the suitability of this software for any purpose. It is *
19 * provided "as is" without express or implied warranty. *
20 **************************************************************************/
21
22
23#include <TROOT.h>
24#include <TH1F.h>
25#include <TH2F.h>
26#include <THnSparse.h>
27#include <TSystem.h>
f0659f11 28#include <TObjArray.h>
d89b8229 29#include <TInterpreter.h>
30#include <TList.h>
31#include <TLorentzVector.h>
32#include <TRefArray.h>
33#include "TDatabasePDG.h"
34
35#include "AliAnalysisTaskJetBackgroundSubtract.h"
36#include "AliAnalysisManager.h"
37#include "AliAODHandler.h"
38#include "AliAODTrack.h"
39#include "AliAODJet.h"
40#include "AliAODEvent.h"
41#include "AliInputEventHandler.h"
42#include "AliAODJetEventBackground.h"
43
44
45ClassImp(AliAnalysisTaskJetBackgroundSubtract)
46
47AliAnalysisTaskJetBackgroundSubtract::~AliAnalysisTaskJetBackgroundSubtract(){
f0659f11 48 //
49 // destructor
50 //
d89b8229 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),
d89b8229 62 fBackgroundBranch(""),
63 fNonStdFile(""),
64 fReplaceString1("B0"),
65 fReplaceString2("B%d"),
66 fSubtraction(k4Area),
67 fKeepJets(kFALSE),
68 fInJetArrayList(0x0),
69 fOutJetArrayList(0x0),
70 fh2CentvsRho(0x0),
71 fh2CentvsSigma(0x0),
72 fh2MultvsRho(0x0),
73 fh2MultvsSigma(0x0),
74 fh2ShiftEta(0x0),
75 fh2ShiftPhi(0x0),
76 fh2ShiftEtaLeading(0x0),
77 fh2ShiftPhiLeading(0x0),
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),
d89b8229 90 fBackgroundBranch(""),
91 fNonStdFile(""),
92 fReplaceString1("B0"),
93 fReplaceString2("B%d"),
94 fSubtraction(k4Area),
95 fKeepJets(kFALSE),
96 fInJetArrayList(0x0),
97 fOutJetArrayList(0x0),
98 fh2CentvsRho(0x0),
99 fh2CentvsSigma(0x0),
100 fh2MultvsRho(0x0),
101 fh2MultvsSigma(0x0),
102 fh2ShiftEta(0x0),
103 fh2ShiftPhi(0x0),
104 fh2ShiftEtaLeading(0x0),
105 fh2ShiftPhiLeading(0x0),
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
d89b8229 117 //
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++){
d89b8229 128 TObjString *ostr = (TObjString*)fJBArray->At(iJB);
129
130
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;
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
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 }
176 return kTRUE;
177}
178
179void AliAnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects()
180{
181
182 //
183 // Create the output container
184 //
185 // Connect the AOD
186
187 if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects() \n");
188
189
190
191 if(fNonStdFile.Length()!=0){
192
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());
196
197 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
198
199 if(!fAODExtension){
200 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
201 }
202 }
203 fAODOut = AODEvent();
204
205 // usually we do not have the input already here
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++){
d89b8229 211 TObjString *ostr = (TObjString*)fJBArray->At(iJB);
212 TString newName(ostr->GetString().Data());
213 if(!newName.Contains(fReplaceString1.Data())){
214 Printf("%s:%d cannot replace string %s in %s",(char*)__FILE__,__LINE__,fReplaceString1.Data(),
215 newName.Data());
216 continue;
217 }
218
219
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);
223 newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
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
229 if(fDebug){
230 Printf("%s:%d created branch \n %s from \n %s",(char*)__FILE__,__LINE__,newName.Data(),
231 ostr->GetString().Data());
232 }
233 tca->SetName(newName.Data());
234 AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
235 }
236
237
238 if(!fHistList)fHistList = new TList();
239 fHistList->SetOwner();
240 PostData(1, fHistList); // post data in any case once
241
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++){
d89b8229 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));
260
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.);
262 fHistList->Add(hTmp);
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);
265 }
266
267 Bool_t oldStatus = TH1::AddDirectoryStatus();
268 TH1::AddDirectory(kFALSE);
269
270 //
271 // Histogram booking, add som control histograms here
272 //
273
274
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.);
277 fHistList->Add(fh2CentvsRho);
278 fHistList->Add(fh2CentvsSigma);
279
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.);
282 fHistList->Add(fh2MultvsRho);
283 fHistList->Add(fh2MultvsSigma);
284
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
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)
d89b8229 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 //
d89b8229 327
328 if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserExec() \n");
329 ResetOutJets();
f0659f11 330 if(fBackgroundBranch.Length()==0||(fJBArray?fJBArray->GetEntries():0)==0){
d89b8229 331 if(fDebug)Printf("%s:%d No background subtraction done",(char*)__FILE__,__LINE__);
332 PostData(1,fHistList);
333 }
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
339
340
341 AliAODJetEventBackground* evBkg = 0;
342 TClonesArray* bkgClusters = 0;
343 TClonesArray* bkgClustersRC = 0;
344 TString bkgClusterName(fBackgroundBranch.Data());
345 bkgClusterName.ReplaceAll(Form("%s_",AliAODJetEventBackground::StdBranchName()),"");
346 TString bkgClusterRCName(Form("%s%s",bkgClusterName.Data(),"RandomCone"));
347
348 if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODOut){
349 evBkg = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
350 bkgClusters = (TClonesArray*)(fAODOut->FindListObject(bkgClusterName.Data()));
351 bkgClustersRC = (TClonesArray*)(fAODOut->FindListObject(bkgClusterRCName.Data()));
352
353 if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
354 if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
355 if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
356 }
357 if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODExtension){
358 evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
359 bkgClusters = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(bkgClusterName.Data()));
360 bkgClustersRC = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(bkgClusterRCName.Data()));
361 if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
362 if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
363
364 if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
365 }
366
367 if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODIn){
368 evBkg = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
369 bkgClusters = (TClonesArray*)(fAODIn->FindListObject(bkgClusterName.Data()));
370 bkgClustersRC = (TClonesArray*)(fAODIn->FindListObject(bkgClusterRCName.Data()));
371
372 if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
373 if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
374 if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
375 }
376
377 if(!evBkg&&(fSubtraction==kArea||fSubtraction==kRhoRecalc||fSubtraction==k4Area)){
378 if(fDebug){
379 Printf("%s:%d Backgroundbranch %s not found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
380 PrintAODContents();
381 }
382 PostData(1,fHistList);
383 return;
384 }
385
386 if(!bkgClusters&&(fSubtraction==kRhoRecalc)){
387 if(fDebug){
388 Printf("%s:%d Background cluster branch %s not found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
389 PrintAODContents();
390 }
391 PostData(1,fHistList);
392 return;
393 }
394
395 if(!bkgClustersRC&&(fSubtraction==kRhoRC)){
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 }
403 // LOOP over all jet branches and subtract the background
404
405 Float_t rho = 0;
406 Float_t sigma=0.;
407 Double_t meanarea = 0;
408 TLorentzVector backgroundv;
409 Float_t cent=0.;
410
411 if(fAODOut)cent = fAODOut->GetHeader()->GetCentrality();
412 if(fAODIn) cent = fAODIn->GetHeader()->GetCentrality();
413
414 if(evBkg)sigma=evBkg->GetSigma(1);
415
416 if(fSubtraction==kArea) rho = evBkg->GetBackground(1);
417 if(fSubtraction==k4Area){
418 rho = evBkg->GetBackground(0);
419 sigma=evBkg->GetSigma(0);
420 }
421
422 if(fSubtraction==kRhoRecalc){
423 meanarea=evBkg->GetMeanarea(1);
424 rho =RecalcRho(bkgClusters,meanarea);
425 }
426 if(fSubtraction==kRhoRC) rho=RhoRC(bkgClustersRC);
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
440 fh2CentvsRho->Fill(cent,rho);
441 fh2CentvsSigma->Fill(cent,sigma);
442
443 fh2MultvsRho->Fill(mult,rho);
444 fh2MultvsSigma->Fill(mult,sigma);
445
446 for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
447 TClonesArray* jarray = (TClonesArray*)fInJetArrayList->At(iJB);
448 TClonesArray* jarrayOut = (TClonesArray*)fOutJetArrayList->At(iJB);
449
450 if(!jarray||!jarrayOut){
451 Printf("%s:%d Array not found %d: %p %p",(char*)__FILE__,__LINE__,iJB,jarray,jarrayOut);
452 continue;
453 }
454 TH2F* h2PtInOut = (TH2F*)fHistList->FindObject(Form("h2PtInPtOut_%d",iJB));
455 THnSparseF* hnDPtAreaCentMult = (THnSparseF*)fHistList->FindObject(Form("hnDPtAreaCentMult_%d",iJB));
456 // loop over all jets
457 Int_t nOut = 0;
458
459 Double_t deltaPt[4];
460 deltaPt[2] = cent;
461 deltaPt[3] = mult;
462
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;
467 Float_t ptSub = 0;
468
469
470 if(fSubtraction==kArea){
471 Double_t background = rho * jet->EffectiveAreaCharged();
472 ptSub = jet->Pt() - background;
473 if(fDebug>2){
474 Printf("%s:%d Jet %d %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub);
475 }
476 if(ptSub<=0){
477 // optionally rescale it and keep??
478 if(fKeepJets){
479 bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
480 }
481 else{
482 bAdd = false;
483 }
484 }
485 else{
486 bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
487 }
488 // add background estimates to the new jet object
489 // allows to recover old p_T and rho...
490 tmpNewJet.SetBgEnergy(background,0);
491 tmpNewJet.SetPtSubtracted(ptSub,0);
492 }// kAREA
493 else if(fSubtraction==kRhoRecalc){
494 Double_t background = rho * jet->EffectiveAreaCharged();
495 ptSub = jet->Pt() - background;
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);}
498 if(ptSub<=0){
499 // optionally rescale it and keep
500 if(fKeepJets){
501 bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
502 }
503 else{
504 bAdd = false;
505 }
506 }
507 else{
508 bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
509 }
510 // add background estimates to the new jet object
511 // allows to recover old p_T and rho...
512 tmpNewJet.SetBgEnergy(background,0);
513 tmpNewJet.SetPtSubtracted(ptSub,0);
514 }//kRhoRecalc
515 else if(fSubtraction==kRhoRC){
516 Double_t background = rho * jet->EffectiveAreaCharged();
517 ptSub = jet->Pt() - background;
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);}
519 if(ptSub<=0){
520 if(fKeepJets){
521 bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
522 }
523 else{
524 bAdd = false;
525 }
526 }
527 else{
528 bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
529 }
530 // add background estimates to the new jet object
531 // allows to recover old p_T and rho...
532 tmpNewJet.SetBgEnergy(background,0);
533 tmpNewJet.SetPtSubtracted(ptSub,0);
534 }//kRhoRC
535
536 else if(fSubtraction==k4Area&&jet->VectorAreaCharged()){
537 backgroundv.SetPxPyPzE(rho*(jet->VectorAreaCharged())->Px(),rho*(jet->VectorAreaCharged())->Py(),rho*(jet->VectorAreaCharged())->Pz(),rho*(jet->VectorAreaCharged())->E());
538 ptSub = jet->Pt()-backgroundv.Pt();
539 if((backgroundv.E()>=jet->E())||(backgroundv.Pt()>=jet->Pt())){
540 if(fKeepJets){
541 bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
542 }
543 else{
544 bAdd = false;
545 }
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...
552 tmpNewJet.SetBgEnergy(backgroundv.Pt(),0);
553 tmpNewJet.SetPtSubtracted(ptSub,0);
554
555 }//kArea4vector
556
557 if(bAdd){
558 AliAODJet *newJet = new ((*jarrayOut)[nOut++]) AliAODJet(tmpNewJet);
559 // what about track references, clear for now...
560 if(fSubtraction==k4Area){
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
566 // set the references
567 newJet->GetRefTracks()->Clear();
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 }
573 }
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 }
580 }
581 if(jarrayOut)jarrayOut->Sort();
582 }
583
584 PostData(1, fHistList);
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));
602
603
604
605 return kTRUE;
606}
607
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
623Double_t AliAnalysisTaskJetBackgroundSubtract::RecalcRho(TClonesArray* bkgClusters,Double_t meanarea){
624
f0659f11 625 //
626 // recalc rhoo
627 //
628
d89b8229 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.
d89b8229 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;
d89b8229 655 ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged();
656 count=count+1;}
657 if(count!=0) rho=ptarea/count;
658 }
659 return rho;
660}
661
f0659f11 662Double_t AliAnalysisTaskJetBackgroundSubtract::RhoRC(TClonesArray* bkgClustersRC){
663
664 //
665 // calc rho from random cones
666 //
d89b8229 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.
d89b8229 672 TClonesArray* jarray=0;
673 for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
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
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;
d89b8229 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
707void AliAnalysisTaskJetBackgroundSubtract::ResetOutJets(){
f0659f11 708 //
709 // Reset the output jets
710 //
711
d89b8229 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}
718
719
720void AliAnalysisTaskJetBackgroundSubtract::PrintAODContents(){
f0659f11 721
722 //
723 // guess from the name what this function does
724 //
725
d89b8229 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}
739
740Int_t AliAnalysisTaskJetBackgroundSubtract::MultFromJetRefs(TClonesArray* jets){
f0659f11 741 //
742 // calculate multiplicty based on jet references
743 //
744
d89b8229 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}