1 // **************************************
2 // Task used for the correction of determiantion of reconstructed jet spectra
3 // Compares input (gen) and output (rec) jets
4 // *******************************************
7 /**************************************************************************
8 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
10 * Author: The ALICE Off-line Project. *
11 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
26 #include <THnSparse.h>
28 #include <TObjArray.h>
29 #include <TInterpreter.h>
31 #include <TLorentzVector.h>
32 #include <TRefArray.h>
33 #include "TDatabasePDG.h"
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"
45 ClassImp(AliAnalysisTaskJetBackgroundSubtract)
47 AliAnalysisTaskJetBackgroundSubtract::~AliAnalysisTaskJetBackgroundSubtract(){
52 delete fOutJetArrayList;
53 delete fInJetArrayList;
56 AliAnalysisTaskJetBackgroundSubtract::AliAnalysisTaskJetBackgroundSubtract():
62 fBackgroundBranch(""),
64 fReplaceString1("B0"),
65 fReplaceString2("B%d"),
69 fOutJetArrayList(0x0),
76 fh2ShiftEtaLeading(0x0),
77 fh2ShiftPhiLeading(0x0),
83 AliAnalysisTaskJetBackgroundSubtract::AliAnalysisTaskJetBackgroundSubtract(const char* name):
85 AliAnalysisTaskSE(name),
90 fBackgroundBranch(""),
92 fReplaceString1("B0"),
93 fReplaceString2("B%d"),
97 fOutJetArrayList(0x0),
104 fh2ShiftEtaLeading(0x0),
105 fh2ShiftPhiLeading(0x0),
108 DefineOutput(1, TList::Class());
113 Bool_t AliAnalysisTaskJetBackgroundSubtract::Notify()
116 // exec with every new file
118 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
122 // Now we also have the Input Event Available! Fillvthe pointers in the list
124 fInJetArrayList->Clear();
125 fOutJetArrayList->Clear();
127 for(int iJB = 0;iJB<(fJBArray?fJBArray->GetEntries():0);iJB++){
128 TObjString *ostr = (TObjString*)fJBArray->At(iJB);
131 TClonesArray* jarray = 0;
132 if(!jarray&&fAODOut){
133 jarray = (TClonesArray*)(fAODOut->FindListObject(ostr->GetString().Data()));
135 if(!jarray&&fAODExtension){
136 jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(ostr->GetString().Data()));
139 jarray = (TClonesArray*)(fAODIn->FindListObject(ostr->GetString().Data()));
144 Printf("%s:%d Input jet branch %s not found",(char*)__FILE__,__LINE__,ostr->GetString().Data());
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());
157 if(!jarrayOut&&fAODOut){
158 jarrayOut = (TClonesArray*)(fAODOut->FindListObject(newName.Data()));
160 if(!jarrayOut&&fAODExtension){
161 jarrayOut = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(newName.Data()));
166 Printf("%s:%d Output jet branch %s not found",(char*)__FILE__,__LINE__,newName.Data());
171 if(jarrayOut&&jarray){
172 fOutJetArrayList->Add(jarrayOut);
173 fInJetArrayList->Add(jarray);
179 void AliAnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects()
183 // Create the output container
187 if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects() \n");
191 if(fNonStdFile.Length()!=0){
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());
197 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
200 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
203 fAODOut = AODEvent();
205 // usually we do not have the input already here
207 if(!fInJetArrayList)fInJetArrayList =new TList();
208 if(!fOutJetArrayList)fOutJetArrayList =new TList();
210 for(int iJB = 0;iJB<(fJBArray?fJBArray->GetEntries():0);iJB++){
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(),
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());
230 Printf("%s:%d created branch \n %s from \n %s",(char*)__FILE__,__LINE__,newName.Data(),
231 ostr->GetString().Data());
233 tca->SetName(newName.Data());
234 AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
238 if(!fHistList)fHistList = new TList();
239 fHistList->SetOwner();
240 PostData(1, fHistList); // post data in any case once
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};
250 for(int iJB = 0;iJB<(fJBArray?fJBArray->GetEntries():0);iJB++){
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(),
259 newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
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);
267 Bool_t oldStatus = TH1::AddDirectoryStatus();
268 TH1::AddDirectory(kFALSE);
271 // Histogram booking, add som control histograms here
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);
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);
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);
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));
303 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
306 TH1::AddDirectory(oldStatus);
308 if(fBackgroundBranch.Length()==0)
309 AliError(Form("%s:%d No BackgroundBranch defined",(char*)__FILE__,__LINE__));
310 if((fJBArray?fJBArray->GetEntries():0)==0)
311 AliError(Form("%s:%d No Jet Branches defined defined",(char*)__FILE__,__LINE__));
314 void AliAnalysisTaskJetBackgroundSubtract::Init()
319 if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::Init() \n");
322 void AliAnalysisTaskJetBackgroundSubtract::UserExec(Option_t */*option*/)
325 // Execute for every selected event
328 if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserExec() \n");
330 if(fBackgroundBranch.Length()==0||(fJBArray?fJBArray->GetEntries():0)==0){
331 if(fDebug)Printf("%s:%d No background subtraction done",(char*)__FILE__,__LINE__);
332 PostData(1,fHistList);
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);
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"));
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()));
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());
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());
364 if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
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()));
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());
377 if(!evBkg&&(fSubtraction==kArea||fSubtraction==kRhoRecalc||fSubtraction==k4Area)){
379 Printf("%s:%d Backgroundbranch %s not found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
382 PostData(1,fHistList);
386 if(!bkgClusters&&(fSubtraction==kRhoRecalc)){
388 Printf("%s:%d Background cluster branch %s not found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
391 PostData(1,fHistList);
395 if(!bkgClustersRC&&(fSubtraction==kRhoRC)){
397 Printf("%s:%d Background cluster RC branch %s not found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
400 PostData(1,fHistList);
403 // LOOP over all jet branches and subtract the background
407 Double_t meanarea = 0;
408 TLorentzVector backgroundv;
411 if(fAODOut)cent = fAODOut->GetHeader()->GetCentrality();
412 if(fAODIn) cent = fAODIn->GetHeader()->GetCentrality();
414 if(evBkg)sigma=evBkg->GetSigma(1);
416 if(fSubtraction==kArea) rho = evBkg->GetBackground(1);
417 if(fSubtraction==k4Area){
418 rho = evBkg->GetBackground(0);
419 sigma=evBkg->GetSigma(0);
422 if(fSubtraction==kRhoRecalc){
423 meanarea=evBkg->GetMeanarea(1);
424 rho =RecalcRho(bkgClusters,meanarea);
426 if(fSubtraction==kRhoRC) rho=RhoRC(bkgClustersRC);
429 for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
430 TClonesArray* jarray = (TClonesArray*)fInJetArrayList->At(iJB);
432 TString tmp(jarray->GetName());
433 if(tmp.Contains("cluster")){
434 mult = MultFromJetRefs(jarray);
440 fh2CentvsRho->Fill(cent,rho);
441 fh2CentvsSigma->Fill(cent,sigma);
443 fh2MultvsRho->Fill(mult,rho);
444 fh2MultvsSigma->Fill(mult,sigma);
446 for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
447 TClonesArray* jarray = (TClonesArray*)fInJetArrayList->At(iJB);
448 TClonesArray* jarrayOut = (TClonesArray*)fOutJetArrayList->At(iJB);
450 if(!jarray||!jarrayOut){
451 Printf("%s:%d Array not found %d: %p %p",(char*)__FILE__,__LINE__,iJB,jarray,jarrayOut);
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
463 for(int i = 0;i < jarray->GetEntriesFast();i++){
464 AliAODJet *jet = (AliAODJet*)jarray->At(i);
465 AliAODJet tmpNewJet(*jet);
470 if(fSubtraction==kArea){
471 Double_t background = rho * jet->EffectiveAreaCharged();
472 ptSub = jet->Pt() - background;
474 Printf("%s:%d Jet %d %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub);
477 // optionally rescale it and keep??
479 bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
486 bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
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);
493 else if(fSubtraction==kRhoRecalc){
494 Double_t background = rho * jet->EffectiveAreaCharged();
495 ptSub = jet->Pt() - background;
497 Printf("%s:%d Jet %d %3.3f %3.3f %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub,background,rho);}
499 // optionally rescale it and keep
501 bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
508 bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
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);
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);}
521 bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
528 bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
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);
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())){
541 bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
548 bAdd = RescaleJet4vector(&tmpNewJet,backgroundv);
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);
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());}}
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);
574 if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
575 if(hnDPtAreaCentMult){
577 deltaPt[1] = jet->EffectiveAreaCharged();
578 hnDPtAreaCentMult->Fill(deltaPt);
581 if(jarrayOut)jarrayOut->Sort();
584 PostData(1, fHistList);
587 void AliAnalysisTaskJetBackgroundSubtract::Terminate(Option_t */*option*/)
589 // Terminate analysis
591 if (fDebug > 1) printf("AnalysisJetBackgroundSubtract: Terminate() \n");
594 Bool_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));
608 Bool_t AliAnalysisTaskJetBackgroundSubtract::RescaleJet4vector(AliAODJet *jet,TLorentzVector backgroundv){
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());
623 Double_t AliAnalysisTaskJetBackgroundSubtract::RecalcRho(TClonesArray* bkgClusters,Double_t meanarea){
632 const Double_t rLimit2=0.8*0.8; //2*jet radius.
633 TClonesArray* jarray=0;
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);
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());
654 if((distance1<rLimit2)||(distance2<rLimit2)) continue;
655 ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged();
657 if(count!=0) rho=ptarea/count;
662 Double_t AliAnalysisTaskJetBackgroundSubtract::RhoRC(TClonesArray* bkgClustersRC){
665 // calc rho from random cones
671 const Double_t rLimit2=0.8*0.8; //2*jet radius.
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);
680 if(!jarray)return rho;
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());
692 if((distance1<rLimit2)||(distance2<rLimit2)) continue;
693 ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged();
695 if(count!=0) rho=ptarea/count; }
707 void AliAnalysisTaskJetBackgroundSubtract::ResetOutJets(){
709 // Reset the output jets
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();
720 void AliAnalysisTaskJetBackgroundSubtract::PrintAODContents(){
723 // guess from the name what this function does
727 Printf("%s:%d >>>>>> Input",(char*)__FILE__,__LINE__);
731 Printf("%s:%d >>>>>> Extenstion",(char*)__FILE__,__LINE__);
732 fAODExtension->GetAOD()->Print();
735 Printf("%s:%d >>>>>> Output",(char*)__FILE__,__LINE__);
740 Int_t AliAnalysisTaskJetBackgroundSubtract::MultFromJetRefs(TClonesArray* jets){
742 // calculate multiplicty based on jet references
748 for(int ij = 0;ij < jets->GetEntries();++ij){
749 AliAODJet* jet = (AliAODJet*)jets->At(ij);
751 TRefArray *refs = jet->GetRefTracks();
753 refMult += refs->GetEntries();
759 void AliAnalysisTaskJetBackgroundSubtract::AddJetBranch(const char* c){
760 if(!fJBArray)fJBArray = new TObjArray();
761 fJBArray->Add(new TObjString(c));