#include <TH2F.h>
#include <THnSparse.h>
#include <TSystem.h>
+#include <TObjArray.h>
#include <TInterpreter.h>
#include <TList.h>
#include <TLorentzVector.h>
ClassImp(AliAnalysisTaskJetBackgroundSubtract)
AliAnalysisTaskJetBackgroundSubtract::~AliAnalysisTaskJetBackgroundSubtract(){
+ //
+ // destructor
+ //
delete fJBArray;
delete fOutJetArrayList;
delete fInJetArrayList;
fAODOut(0x0),
fAODIn(0x0),
fAODExtension(0x0),
- fJBArray(new TObjArray()),
+ fJBArray(0),
fBackgroundBranch(""),
fNonStdFile(""),
fReplaceString1("B0"),
fReplaceString2("B%d"),
fSubtraction(k4Area),
+ fKeepJets(kFALSE),
fInJetArrayList(0x0),
fOutJetArrayList(0x0),
fh2CentvsRho(0x0),
fh2CentvsSigma(0x0),
+ fh2MultvsRho(0x0),
+ fh2MultvsSigma(0x0),
fh2ShiftEta(0x0),
fh2ShiftPhi(0x0),
fh2ShiftEtaLeading(0x0),
fAODOut(0x0),
fAODIn(0x0),
fAODExtension(0x0),
- fJBArray(new TObjArray()),
+ fJBArray(0),
fBackgroundBranch(""),
fNonStdFile(""),
fReplaceString1("B0"),
fReplaceString2("B%d"),
fSubtraction(k4Area),
+ fKeepJets(kFALSE),
fInJetArrayList(0x0),
fOutJetArrayList(0x0),
fh2CentvsRho(0x0),
fh2CentvsSigma(0x0),
+ fh2MultvsRho(0x0),
+ fh2MultvsSigma(0x0),
fh2ShiftEta(0x0),
fh2ShiftPhi(0x0),
fh2ShiftEtaLeading(0x0),
Bool_t AliAnalysisTaskJetBackgroundSubtract::Notify()
{
+ //
+ // exec with every new file
//
fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
fInJetArrayList->Clear();
fOutJetArrayList->Clear();
- for(int iJB = 0;iJB<fJBArray->GetEntries();iJB++){
+ for(int iJB = 0;iJB<(fJBArray?fJBArray->GetEntries():0);iJB++){
TObjString *ostr = (TObjString*)fJBArray->At(iJB);
TString newName(ostr->GetString().Data());
newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
TClonesArray* jarrayOut = 0;
+ if(newName.CompareTo(ostr->GetString())==0){
+ Printf("%s:%d Input and output branch would have the same name, skipping %s ",(char*)__FILE__,__LINE__,ostr->GetString().Data());
+ continue;
+ }
+
if(!jarrayOut&&fAODOut){
jarrayOut = (TClonesArray*)(fAODOut->FindListObject(newName.Data()));
}
if(!fInJetArrayList)fInJetArrayList =new TList();
if(!fOutJetArrayList)fOutJetArrayList =new TList();
- for(int iJB = 0;iJB<fJBArray->GetEntries();iJB++){
+ for(int iJB = 0;iJB<(fJBArray?fJBArray->GetEntries():0);iJB++){
TObjString *ostr = (TObjString*)fJBArray->At(iJB);
TString newName(ostr->GetString().Data());
if(!newName.Contains(fReplaceString1.Data())){
continue;
}
+
// add a new branch to the output for the background subtracted jets take the names from
// the input jets and replace the background flag names
TClonesArray *tca = new TClonesArray("AliAODJet", 0);
newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
+ if(newName.CompareTo(ostr->GetString())==0){
+ Printf("%s:%d Input and output branch would have the same name, skipping: %s ",(char*)__FILE__,__LINE__,ostr->GetString().Data());
+ continue;
+ }
+
if(fDebug){
Printf("%s:%d created branch \n %s from \n %s",(char*)__FILE__,__LINE__,newName.Data(),
ostr->GetString().Data());
fHistList->SetOwner();
PostData(1, fHistList); // post data in any case once
- for(int iJB = 0;iJB<fJBArray->GetEntries();iJB++){
+ //
+
+ // delta pT vs. area vs. cent vs. mult
+ const Int_t nSparseBinsDelta = 4;
+ const Int_t nBinsDelta[nSparseBinsDelta] = { 241, 10, 10, 25};
+ const Double_t xminDelta[nSparseBinsDelta] = {-120.5, 0, 0, 0};
+ const Double_t xmaxDelta[nSparseBinsDelta] = { 120.5, 1.0, 100,5000};
+
+ for(int iJB = 0;iJB<(fJBArray?fJBArray->GetEntries():0);iJB++){
TObjString *ostr = (TObjString*)fJBArray->At(iJB);
TString oldName(ostr->GetString().Data());
TString newName(ostr->GetString().Data());
continue;
}
newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
+
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.);
fHistList->Add(hTmp);
+ THnSparseF *hFTmp = new THnSparseF(Form("hnDPtAreaCentMult_%d",iJB),Form("%s delta;#delta p_{T};Area;cent;mult",newName.Data()),nSparseBinsDelta,nBinsDelta,xminDelta,xmaxDelta);
+ fHistList->Add(hFTmp);
}
Bool_t oldStatus = TH1::AddDirectoryStatus();
//
- fh2CentvsRho = new TH2F("fh2CentvsRho","centrality vs background density",100,0.,100.,100,0.,150.);
- fh2CentvsSigma = new TH2F("fh2CentvsSigma","centrality vs backgroun sigma",100,0.,100.,100,0.,150.);
+ fh2CentvsRho = new TH2F("fh2CentvsRho","centrality vs background density", 100,0.,100.,600,0.,300.);
+ fh2CentvsSigma = new TH2F("fh2CentvsSigma","centrality vs backgroun sigma",100,0.,100.,500,0.,50.);
fHistList->Add(fh2CentvsRho);
fHistList->Add(fh2CentvsSigma);
+ fh2MultvsRho = new TH2F("fh2MultvsRho","mult vs background density", 100,0.,5000.,600,0.,300.);
+ fh2MultvsSigma = new TH2F("fh2MultvsSigma","mult vs background sigma",100,0.,5000.,500,0.,50.);
+ fHistList->Add(fh2MultvsRho);
+ fHistList->Add(fh2MultvsSigma);
+
if(fSubtraction==k4Area){
fh2ShiftEta = new TH2F("fh2ShiftEta","extended correction Eta",100,-0.9,0.9,100,-0.9,0.9);
fh2ShiftPhi = new TH2F("fh2ShiftPhi","extended correction Phi",100,0.,6.5,100,0.,6.5);
fHistList->Add(fh2ShiftPhiLeading);
}
-
// =========== Switch on Sumw2 for all histos ===========
for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
if(fBackgroundBranch.Length()==0)
AliError(Form("%s:%d No BackgroundBranch defined",(char*)__FILE__,__LINE__));
- if(fJBArray->GetEntries()==0)
+ if((fJBArray?fJBArray->GetEntries():0)==0)
AliError(Form("%s:%d No Jet Branches defined defined",(char*)__FILE__,__LINE__));
}
void AliAnalysisTaskJetBackgroundSubtract::UserExec(Option_t */*option*/)
{
+ //
+ // Execute for every selected event
+ //
if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserExec() \n");
ResetOutJets();
- if(fBackgroundBranch.Length()==0||fJBArray->GetEntries()==0){
+ if(fBackgroundBranch.Length()==0||(fJBArray?fJBArray->GetEntries():0)==0){
if(fDebug)Printf("%s:%d No background subtraction done",(char*)__FILE__,__LINE__);
PostData(1,fHistList);
}
Double_t meanarea = 0;
TLorentzVector backgroundv;
Float_t cent=0.;
+
if(fAODOut)cent = fAODOut->GetHeader()->GetCentrality();
+ if(fAODIn) cent = fAODIn->GetHeader()->GetCentrality();
+
if(evBkg)sigma=evBkg->GetSigma(1);
if(fSubtraction==kArea) rho = evBkg->GetBackground(1);
rho =RecalcRho(bkgClusters,meanarea);
}
if(fSubtraction==kRhoRC) rho=RhoRC(bkgClustersRC);
+
+ Float_t mult = 0;
+ for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
+ TClonesArray* jarray = (TClonesArray*)fInJetArrayList->At(iJB);
+ if(jarray){
+ TString tmp(jarray->GetName());
+ if(tmp.Contains("cluster")){
+ mult = MultFromJetRefs(jarray);
+ if(mult>0)break;
+ }
+ }
+ }
+
fh2CentvsRho->Fill(cent,rho);
fh2CentvsSigma->Fill(cent,sigma);
+
+ fh2MultvsRho->Fill(mult,rho);
+ fh2MultvsSigma->Fill(mult,sigma);
for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
TClonesArray* jarray = (TClonesArray*)fInJetArrayList->At(iJB);
continue;
}
TH2F* h2PtInOut = (TH2F*)fHistList->FindObject(Form("h2PtInPtOut_%d",iJB));
+ THnSparseF* hnDPtAreaCentMult = (THnSparseF*)fHistList->FindObject(Form("hnDPtAreaCentMult_%d",iJB));
// loop over all jets
Int_t nOut = 0;
+ Double_t deltaPt[4];
+ deltaPt[2] = cent;
+ deltaPt[3] = mult;
for(int i = 0;i < jarray->GetEntriesFast();i++){
AliAODJet *jet = (AliAODJet*)jarray->At(i);
AliAODJet tmpNewJet(*jet);
Bool_t bAdd = false;
-
+ Float_t ptSub = 0;
+
if(fSubtraction==kArea){
Double_t background = rho * jet->EffectiveAreaCharged();
- Float_t ptSub = jet->Pt() - background;
+ ptSub = jet->Pt() - background;
if(fDebug>2){
Printf("%s:%d Jet %d %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub);
}
- if(ptSub<0){
+ if(ptSub<=0){
// optionally rescale it and keep??
- bAdd = false; // RescaleJetMomentum(&tmpNewJet,0.1);
- if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
+ if(fKeepJets){
+ bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
+ }
+ else{
+ bAdd = false;
+ }
}
else{
bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
- if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
}
// add background estimates to the new jet object
// allows to recover old p_T and rho...
tmpNewJet.SetBgEnergy(background,0);
+ tmpNewJet.SetPtSubtracted(ptSub,0);
}// kAREA
else if(fSubtraction==kRhoRecalc){
Double_t background = rho * jet->EffectiveAreaCharged();
- Float_t ptSub = jet->Pt() - background;
+ ptSub = jet->Pt() - background;
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);}
- if(ptSub<0){
- // optionally rescale it and keep??
- bAdd = false;// RescaleJetMomentum(&tmpNewJet,0.1);
- if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
+ if(ptSub<=0){
+ // optionally rescale it and keep
+ if(fKeepJets){
+ bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
+ }
+ else{
+ bAdd = false;
+ }
}
else{
bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
- if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
}
// add background estimates to the new jet object
// allows to recover old p_T and rho...
tmpNewJet.SetBgEnergy(background,0);
-
+ tmpNewJet.SetPtSubtracted(ptSub,0);
}//kRhoRecalc
else if(fSubtraction==kRhoRC){
Double_t background = rho * jet->EffectiveAreaCharged();
- Float_t ptSub = jet->Pt() - background;
+ ptSub = jet->Pt() - background;
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);}
- if(ptSub<0){
- // optionally rescale it and keep??
- bAdd = false; // RescaleJetMomentum(&tmpNewJet,0.1);
- if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
+ if(ptSub<=0){
+ if(fKeepJets){
+ bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
+ }
+ else{
+ bAdd = false;
+ }
}
else{
bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
- if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
}
// add background estimates to the new jet object
// allows to recover old p_T and rho...
tmpNewJet.SetBgEnergy(background,0);
-
+ tmpNewJet.SetPtSubtracted(ptSub,0);
}//kRhoRC
else if(fSubtraction==k4Area&&jet->VectorAreaCharged()){
-
-
backgroundv.SetPxPyPzE(rho*(jet->VectorAreaCharged())->Px(),rho*(jet->VectorAreaCharged())->Py(),rho*(jet->VectorAreaCharged())->Pz(),rho*(jet->VectorAreaCharged())->E());
- if((backgroundv.E()>jet->E())&&(backgroundv.Pt()>jet->Pt())){
-
- // optionally rescale it and keep??
- bAdd = false; // RescaleJetMomentum(&tmpNewJet,0.1);
- if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),jet->Pt()-backgroundv.Pt());
+ ptSub = jet->Pt()-backgroundv.Pt();
+ if((backgroundv.E()>=jet->E())||(backgroundv.Pt()>=jet->Pt())){
+ if(fKeepJets){
+ bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
+ }
+ else{
+ bAdd = false;
+ }
}
else{
bAdd = RescaleJet4vector(&tmpNewJet,backgroundv);
}
// add background estimates to the new jet object
// allows to recover old p_T and rho...
- tmpNewJet.SetBgEnergy(backgroundv.P(),0);
+ tmpNewJet.SetBgEnergy(backgroundv.Pt(),0);
+ tmpNewJet.SetPtSubtracted(ptSub,0);
}//kArea4vector
-
-
-
-
-
if(bAdd){
AliAODJet *newJet = new ((*jarrayOut)[nOut++]) AliAODJet(tmpNewJet);
// what about track references, clear for now...
if(fSubtraction==k4Area){
- if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),jet->Pt()-newJet->Pt());
fh2ShiftEta->Fill(jet->Eta(),newJet->Eta());
fh2ShiftPhi->Fill(jet->Phi(),newJet->Phi());
if(i==0){fh2ShiftEtaLeading->Fill(jet->Eta(),newJet->Eta());
if(vp)newJet->AddTrack(vp);
}
}
+ if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
+ if(hnDPtAreaCentMult){
+ deltaPt[0] = ptSub;
+ deltaPt[1] = jet->EffectiveAreaCharged();
+ hnDPtAreaCentMult->Fill(deltaPt);
+ }
}
- }
- PostData(1, fHistList);
+ if(jarrayOut)jarrayOut->Sort();
+ }
+
+ PostData(1, fHistList);
}
void AliAnalysisTaskJetBackgroundSubtract::Terminate(Option_t */*option*/)
Double_t AliAnalysisTaskJetBackgroundSubtract::RecalcRho(TClonesArray* bkgClusters,Double_t meanarea){
+ //
+ // recalc rhoo
+ //
+
Double_t ptarea=0.;
Int_t count=0;
Double_t rho=0.;
- const Double_t Rlimit2=0.8*0.8; //2*jet radius.
+ const Double_t rLimit2=0.8*0.8; //2*jet radius.
TClonesArray* jarray=0;
for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
(first->Phi()-clus->Phi())*(first->Phi()-clus->Phi());
Double_t distance2= (second->Eta()-clus->Eta())*(second->Eta()-clus->Eta())+
(second->Phi()-clus->Phi())*(second->Phi()-clus->Phi());
- if((distance1<Rlimit2)||(distance2<Rlimit2)) continue;
+ if((distance1<rLimit2)||(distance2<rLimit2)) continue;
ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged();
count=count+1;}
if(count!=0) rho=ptarea/count;
return rho;
}
- Double_t AliAnalysisTaskJetBackgroundSubtract::RhoRC(TClonesArray* bkgClustersRC){
+Double_t AliAnalysisTaskJetBackgroundSubtract::RhoRC(TClonesArray* bkgClustersRC){
+
+ //
+ // calc rho from random cones
+ //
Double_t ptarea=0.;
Int_t count=0;
Double_t rho=0.;
- const Double_t Rlimit2=0.8*0.8; //2*jet radius.
+ const Double_t rLimit2=0.8*0.8; //2*jet radius.
TClonesArray* jarray=0;
for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
TObjString *ostr = (TObjString*)fInJetArrayList->At(iJB);
(first->Phi()-clus->Phi())*(first->Phi()-clus->Phi());
Double_t distance2= (second->Eta()-clus->Eta())*(second->Eta()-clus->Eta())+
(second->Phi()-clus->Phi())*(second->Phi()-clus->Phi());
- if((distance1<Rlimit2)||(distance2<Rlimit2)) continue;
+ if((distance1<rLimit2)||(distance2<rLimit2)) continue;
ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged();
count=count+1;}
if(count!=0) rho=ptarea/count; }
void AliAnalysisTaskJetBackgroundSubtract::ResetOutJets(){
+ //
+ // Reset the output jets
+ //
+
if(!fOutJetArrayList)return;
for(int iJB = 0;iJB<fOutJetArrayList->GetEntries();iJB++){
TClonesArray* jarray = (TClonesArray*)fOutJetArrayList->At(iJB);
void AliAnalysisTaskJetBackgroundSubtract::PrintAODContents(){
+
+ //
+ // guess from the name what this function does
+ //
+
if(fAODIn){
Printf("%s:%d >>>>>> Input",(char*)__FILE__,__LINE__);
fAODIn->Print();
fAODOut->Print();
}
}
+
+Int_t AliAnalysisTaskJetBackgroundSubtract::MultFromJetRefs(TClonesArray* jets){
+ //
+ // calculate multiplicty based on jet references
+ //
+
+ if(!jets)return 0;
+
+ Int_t refMult = 0;
+ for(int ij = 0;ij < jets->GetEntries();++ij){
+ AliAODJet* jet = (AliAODJet*)jets->At(ij);
+ if(!jet)continue;
+ TRefArray *refs = jet->GetRefTracks();
+ if(!refs)continue;
+ refMult += refs->GetEntries();
+ }
+ return refMult;
+
+}
+
+void AliAnalysisTaskJetBackgroundSubtract::AddJetBranch(const char* c){
+ if(!fJBArray)fJBArray = new TObjArray();
+ fJBArray->Add(new TObjString(c));
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff