Added 4 Area backgroudn subtraction

diff --git a/JETAN/AliAnalysisTaskJetBackgroundSubtract.cxx b/JETAN/AliAnalysisTaskJetBackgroundSubtract.cxx
index 5318f5c3aa462ea95fdcc1b9390d2d97d0ba732a..ed0868d39602504898b4894270d76e700eeb96eb 100644 (file)
--- a/JETAN/AliAnalysisTaskJetBackgroundSubtract.cxx
+++ b/JETAN/AliAnalysisTaskJetBackgroundSubtract.cxx
@@ -59,9 +59,15 @@ AliAnalysisTaskJetBackgroundSubtract::AliAnalysisTaskJetBackgroundSubtract():
   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)  
 {
 
@@ -78,9 +84,15 @@ AliAnalysisTaskJetBackgroundSubtract::AliAnalysisTaskJetBackgroundSubtract(const
   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());  
@@ -220,6 +232,24 @@ void AliAnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects()
   //  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));
@@ -327,15 +357,27 @@ void AliAnalysisTaskJetBackgroundSubtract::UserExec(Option_t */*option*/)
   // 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);
     
@@ -352,7 +394,7 @@ void AliAnalysisTaskJetBackgroundSubtract::UserExec(Option_t */*option*/)
       AliAODJet *jet = (AliAODJet*)jarray->At(i);
       AliAODJet tmpNewJet(*jet);
       Bool_t bAdd = false;
-
+     
 
       if(fSubtraction==kArea){ 
        Double_t background = rho * jet->EffectiveAreaCharged();
@@ -411,9 +453,40 @@ void AliAnalysisTaskJetBackgroundSubtract::UserExec(Option_t */*option*/)
 
        }//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();
       }
 
@@ -452,9 +525,27 @@ Bool_t AliAnalysisTaskJetBackgroundSubtract::RescaleJetMomentum(AliAODJet *jet,F
   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.;

diff --git a/JETAN/AliAnalysisTaskJetBackgroundSubtract.h b/JETAN/AliAnalysisTaskJetBackgroundSubtract.h
index 212cd3c9c6a3166434c45c24ac0d717af998dedb..961ced75491ed5e7bf284265124cbecbbbab9afd 100644 (file)
--- a/JETAN/AliAnalysisTaskJetBackgroundSubtract.h
+++ b/JETAN/AliAnalysisTaskJetBackgroundSubtract.h
@@ -66,6 +66,8 @@ class AliAnalysisTaskJetBackgroundSubtract : public AliAnalysisTaskSE
     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();
@@ -83,9 +85,18 @@ class AliAnalysisTaskJetBackgroundSubtract : public AliAnalysisTaskSE
     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

diff --git a/JETAN/AliAnalysisTaskJetCluster.cxx b/JETAN/AliAnalysisTaskJetCluster.cxx
index fc6112685b02cf57708ba2257a6ab2cd25833d2a..a0144c069e84e1750e04169de98750b38a8ff579 100644 (file)
--- a/JETAN/AliAnalysisTaskJetCluster.cxx
+++ b/JETAN/AliAnalysisTaskJetCluster.cxx
@@ -877,7 +877,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
     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
@@ -902,7 +902,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
     }  
     
    
-    
+    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;
@@ -913,6 +913,11 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
        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);
+
+
       }
 
 
@@ -1108,7 +1113,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
      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));    
@@ -1195,11 +1200,14 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
 
  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){
@@ -1249,7 +1257,12 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
        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);
+
+     }