fNonStdFile(""),
fReplaceString1("B0"),
fReplaceString2("B%d"),
- fSubtraction(kRhoRecalc),
+ fSubtraction(k4Area),
fInJetArrayList(0x0),
fOutJetArrayList(0x0),
+ fh2CentvsRho(0x0),
+ fh2CentvsSigma(0x0),
+ fh2ShiftEta(0x0),
+ fh2ShiftPhi(0x0),
+ fh2ShiftEtaLeading(0x0),
+ fh2ShiftPhiLeading(0x0),
fHistList(0x0)
{
fNonStdFile(""),
fReplaceString1("B0"),
fReplaceString2("B%d"),
- fSubtraction(kRhoRecalc),
+ fSubtraction(k4Area),
fInJetArrayList(0x0),
fOutJetArrayList(0x0),
+ fh2CentvsRho(0x0),
+ fh2CentvsSigma(0x0),
+ fh2ShiftEta(0x0),
+ fh2ShiftPhi(0x0),
+ fh2ShiftEtaLeading(0x0),
+ fh2ShiftPhiLeading(0x0),
fHistList(0x0)
{
DefineOutput(1, TList::Class());
// Histogram booking, add som control histograms here
//
+
+ 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.);
+ fHistList->Add(fh2CentvsRho);
+ fHistList->Add(fh2CentvsSigma);
+
+ 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);
+ fh2ShiftEtaLeading = new TH2F("fh2ShiftEtaLeading","extended correction Eta",100,-0.9,0.9,100,-0.9,0.9);
+ fh2ShiftPhiLeading = new TH2F("fh2ShiftPhiLeading","extended correction Phi",100,0.,6.5,100,0.,6.5);
+ fHistList->Add(fh2ShiftEta);
+ fHistList->Add(fh2ShiftPhi);
+ fHistList->Add(fh2ShiftEtaLeading);
+ 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));
// LOOP over all jet branches and subtract the background
Float_t rho = 0;
+ Float_t sigma=0.;
Double_t meanarea = 0;
- if(fSubtraction==kArea)rho = evBkg->GetBackground(2);
+ TLorentzVector backgroundv;
+ Float_t cent=0.;
+ cent = fAODOut->GetHeader()->GetCentrality();
+ sigma=evBkg->GetSigma(1);
+
+ if(fSubtraction==kArea) rho = evBkg->GetBackground(1);
+ if(fSubtraction==k4Area){
+ rho = evBkg->GetBackground(0);
+ sigma=evBkg->GetSigma(0);}
+
if(fSubtraction==kRhoRecalc){
- meanarea=evBkg->GetMeanarea(2);
+ meanarea=evBkg->GetMeanarea(1);
rho =RecalcRho(bkgClusters,meanarea);
- }
- if(fSubtraction==kRhoRC)rho =RhoRC(bkgClustersRC);
+ }
+ if(fSubtraction==kRhoRC) rho=RhoRC(bkgClustersRC);
+ fh2CentvsRho->Fill(cent,rho);
+ fh2CentvsSigma->Fill(cent,sigma);
- for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
+ for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
TClonesArray* jarray = (TClonesArray*)fInJetArrayList->At(iJB);
TClonesArray* jarrayOut = (TClonesArray*)fOutJetArrayList->At(iJB);
AliAODJet *jet = (AliAODJet*)jarray->At(i);
AliAODJet tmpNewJet(*jet);
Bool_t bAdd = false;
-
+
if(fSubtraction==kArea){
Double_t background = rho * jet->EffectiveAreaCharged();
}//kRhoRC
+ else if(fSubtraction==k4Area){
+
+ 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 = RescaleJetMomentum(&tmpNewJet,0.1);
+ if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),0.1);
+ }
+ 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);
+
+ }//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());
+ fh2ShiftPhiLeading->Fill(jet->Phi(),newJet->Phi());}}
+
+
newJet->GetRefTracks()->Clear();
}
Double_t mass = jet->M();
Double_t pNew = jet->P() * scale;
jet->SetPxPyPzE(scale*jet->Px(),scale*jet->Py(),scale*jet->Pz(),TMath::Sqrt(mass*mass+pNew*pNew));
+
+
+
return kTRUE;
}
+Bool_t AliAnalysisTaskJetBackgroundSubtract::RescaleJet4vector(AliAODJet *jet,TLorentzVector backgroundv){
+
+ if(backgroundv.Pt()<0.) return kFALSE;
+ jet->SetPxPyPzE(jet->Px()-backgroundv.Px(),jet->Py()-backgroundv.Py(),jet->Pz()-backgroundv.Pz(),jet->E()-backgroundv.E());
+
+ return kTRUE;
+}
+
+
+
+
+
+
+
+
Double_t AliAnalysisTaskJetBackgroundSubtract::RecalcRho(TClonesArray* bkgClusters,Double_t meanarea){
Double_t ptarea=0.;
AliAnalysisTaskJetBackgroundSubtract(const AliAnalysisTaskJetBackgroundSubtract&);
AliAnalysisTaskJetBackgroundSubtract& operator=(const AliAnalysisTaskJetBackgroundSubtract&);
Bool_t RescaleJetMomentum(AliAODJet *jet,Float_t pT);
+ Bool_t RescaleJet4vector(AliAODJet *jet,TLorentzVector backgroundv);
+
Double_t RecalcRho(TClonesArray* fbkgclusters,Double_t meanarea);
Double_t RhoRC(TClonesArray* fbkgclustersRC);
void ResetOutJets();
Int_t fSubtraction; // Parameter for subtraction mode
TList *fInJetArrayList; //! transient list to make ease the handling of input jets
TList *fOutJetArrayList; //! transient list to make ease the reset of output jets
+
+ TH2F* fh2CentvsRho; //! centrality vs background density
+ TH2F* fh2CentvsSigma; //! centrality vs background sigma
+ TH2F* fh2ShiftEta; //! extended correction Eta
+ TH2F* fh2ShiftPhi; //! extended correction Phi
+ TH2F* fh2ShiftEtaLeading; //! extended correction Eta leading jet
+ TH2F* fh2ShiftPhiLeading; //! extended correction Phi leading jet
+
+
TList *fHistList; //! the histograms output list
- ClassDef(AliAnalysisTaskJetBackgroundSubtract, 3)
+ ClassDef(AliAnalysisTaskJetBackgroundSubtract, 4)
};
#endif
if(externalBackground){
// carefull has to be filled in a task before
// todo, ReArrange to the botom
- pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
+ pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged();
}
pt = leadingJet.Pt() - pTback;
// correlation of leading jet with tracks
}
-
+ TLorentzVector vecareab;
for(int j = 0; j < nRec;j++){
AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
aodOutJet = 0;
aodOutJet = new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
Double_t area1 = clustSeq.area(sortedJets[j]);
aodOutJet->SetEffArea(area1,0);
+ fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
+ vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
+ aodOutJet->SetVectorAreaCharged(&vecareab);
+
+
}
Double_t sigma2=0.;
Double_t meanarea2=0.;
- clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, true);
+ clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
// fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
Int_t nRecOverRan = inclusiveJetsRan.size();
Int_t nRecRan = inclusiveJetsRan.size();
+
if(inclusiveJetsRan.size()>0){
AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
Float_t pt = leadingJet.Pt();
Int_t iCount = 0;
+ TLorentzVector vecarearanb;
+
for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
while(pt<ptCut&&iCount<nRecRan){
aodOutJetRan = new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec);
Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
- aodOutJetRan->SetEffArea(arearan,0); }
+ aodOutJetRan->SetEffArea(arearan,0);
+ fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
+ vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
+ aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
+
+ }