handle reference tracks correctly for background subtraction and randomized events
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskJetBackgroundSubtract.cxx
index ed0868d..a280562 100644 (file)
@@ -168,12 +168,16 @@ void AliAnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects()
   // Connect the AOD
 
   if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects() \n");
+
+
+
   if(fNonStdFile.Length()!=0){
     
     // case that we have an AOD extension we need to fetch the jets from the extended output
     // we identifay the extension aod event by looking for the branchname
     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
-    fAODExtension = aodH->GetExtension(fNonStdFile.Data());
+
+    fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
     
     if(!fAODExtension){
       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
@@ -210,6 +214,7 @@ void AliAnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects()
 
   if(!fHistList)fHistList = new TList();
   fHistList->SetOwner();
+  PostData(1, fHistList); // post data in any case once
 
   for(int iJB = 0;iJB<fJBArray->GetEntries();iJB++){
     TObjString *ostr = (TObjString*)fJBArray->At(iJB);
@@ -221,7 +226,7 @@ void AliAnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects()
       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.,200,0.,200.);
+    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);
   }
 
@@ -328,7 +333,7 @@ void AliAnalysisTaskJetBackgroundSubtract::UserExec(Option_t */*option*/)
     if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
   }
 
-  if(!evBkg&&(fSubtraction==kArea||fSubtraction==kRhoRecalc)){
+  if(!evBkg&&(fSubtraction==kArea||fSubtraction==kRhoRecalc||fSubtraction==k4Area)){
     if(fDebug){
       Printf("%s:%d Backgroundbranch %s not found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
       PrintAODContents();
@@ -361,18 +366,19 @@ void AliAnalysisTaskJetBackgroundSubtract::UserExec(Option_t */*option*/)
    Double_t meanarea = 0;
    TLorentzVector backgroundv;
    Float_t cent=0.;
-   cent = fAODOut->GetHeader()->GetCentrality();
-   sigma=evBkg->GetSigma(1); 
+   if(fAODOut)cent = fAODOut->GetHeader()->GetCentrality();
+   if(evBkg)sigma=evBkg->GetSigma(1); 
 
    if(fSubtraction==kArea) rho = evBkg->GetBackground(1);
    if(fSubtraction==k4Area){
-                rho = evBkg->GetBackground(0);
-                sigma=evBkg->GetSigma(0);}
+     rho = evBkg->GetBackground(0);
+     sigma=evBkg->GetSigma(0);
+   }
 
    if(fSubtraction==kRhoRecalc){
      meanarea=evBkg->GetMeanarea(1);
      rho =RecalcRho(bkgClusters,meanarea);
-        }
+   }
    if(fSubtraction==kRhoRC) rho=RhoRC(bkgClustersRC);
    fh2CentvsRho->Fill(cent,rho);
    fh2CentvsSigma->Fill(cent,sigma);
@@ -382,9 +388,9 @@ void AliAnalysisTaskJetBackgroundSubtract::UserExec(Option_t */*option*/)
     TClonesArray* jarrayOut = (TClonesArray*)fOutJetArrayList->At(iJB);
     
     if(!jarray||!jarrayOut){
-    Printf("%s:%d Array not found %d: %p %p",(char*)__FILE__,__LINE__,iJB,jarray->GetName(),jarrayOut->GetName());
-    continue;
-      }
+      Printf("%s:%d Array not found %d: %p %p",(char*)__FILE__,__LINE__,iJB,jarray,jarrayOut);
+      continue;
+    }
     TH2F* h2PtInOut = (TH2F*)fHistList->FindObject(Form("h2PtInPtOut_%d",iJB));
     // loop over all jets
     Int_t nOut = 0;
@@ -404,8 +410,8 @@ void AliAnalysisTaskJetBackgroundSubtract::UserExec(Option_t */*option*/)
        }
        if(ptSub<0){
          // optionally rescale it and keep??
-         bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
-         if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),0.1);
+         bAdd = false; // RescaleJetMomentum(&tmpNewJet,0.1);
+         if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
        }
        else{
          bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
@@ -422,8 +428,8 @@ void AliAnalysisTaskJetBackgroundSubtract::UserExec(Option_t */*option*/)
          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 = RescaleJetMomentum(&tmpNewJet,0.1);
-         if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),0.1);
+         bAdd = false;// RescaleJetMomentum(&tmpNewJet,0.1);
+         if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
        }
        else{
          bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
@@ -440,8 +446,8 @@ void AliAnalysisTaskJetBackgroundSubtract::UserExec(Option_t */*option*/)
        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 = RescaleJetMomentum(&tmpNewJet,0.1);
-         if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),0.1);
+         bAdd = false; // RescaleJetMomentum(&tmpNewJet,0.1);
+         if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
        }
        else{
          bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
@@ -453,22 +459,23 @@ 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);
+       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());
+        }
+        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  
 
 
@@ -486,26 +493,15 @@ void AliAnalysisTaskJetBackgroundSubtract::UserExec(Option_t */*option*/)
          if(i==0){fh2ShiftEtaLeading->Fill(jet->Eta(),newJet->Eta());
           fh2ShiftPhiLeading->Fill(jet->Phi(),newJet->Phi());}}
 
-
+       // set the references 
        newJet->GetRefTracks()->Clear();
+       TRefArray *refs = jet->GetRefTracks();
+       for(Int_t ir=0;ir<refs->GetEntriesFast();ir++){
+         AliVParticle *vp = dynamic_cast<AliVParticle*>(refs->At(ir));
+         if(vp)newJet->AddTrack(vp);
+       }
       }
-
-
-
-
-
-
     }
-
-
-
-    // subtract the background
-    
-
-    // remove jets??
-
-    // sort jets...
-
   }
   PostData(1, fHistList);
 }
@@ -548,34 +544,37 @@ Bool_t AliAnalysisTaskJetBackgroundSubtract::RescaleJet4vector(AliAODJet *jet,TL
 
 Double_t AliAnalysisTaskJetBackgroundSubtract::RecalcRho(TClonesArray* bkgClusters,Double_t meanarea){
   
-       Double_t ptarea=0.;
-       Int_t count=0;
-       Double_t rho=0.; 
-       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);
-        TString jetref=ostr->GetString().Data();
-        if(jetref.Contains("ANTIKT04")){ 
-         jarray = (TClonesArray*)fInJetArrayList->At(iJB);}}
-       if(jarray->GetEntries()>=2){ 
-         AliAODJet *first = (AliAODJet*)(jarray->At(0)); 
-         AliAODJet *second= (AliAODJet*)(jarray->At(1)); 
-         for(Int_t k=0;k<bkgClusters->GetEntriesFast();k++){
-           AliAODJet *clus = (AliAODJet*)(bkgClusters->At(k));
-           if(TMath::Abs(clus->Eta())>0.5) continue;
-           if((clus->EffectiveAreaCharged())<0.1*meanarea) continue; 
-           Double_t distance1=(first->Eta()-clus->Eta())*(first->Eta()-clus->Eta())+
-             (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;    
-           ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged(); 
-           count=count+1;}
-         if(count!=0) rho=ptarea/count; 
-       }        
-       return rho;
+  Double_t ptarea=0.;
+  Int_t count=0;
+  Double_t rho=0.; 
+  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);
+    TString jetref=ostr->GetString().Data();
+    if(jetref.Contains("ANTIKT04")){ 
+      jarray = (TClonesArray*)fInJetArrayList->At(iJB);
+    }
+  }
+  if(!jarray)return rho;
+  if(jarray->GetEntries()>=2){ 
+    AliAODJet *first = (AliAODJet*)(jarray->At(0)); 
+    AliAODJet *second= (AliAODJet*)(jarray->At(1)); 
+    for(Int_t k=0;k<bkgClusters->GetEntriesFast();k++){
+      AliAODJet *clus = (AliAODJet*)(bkgClusters->At(k));
+      if(TMath::Abs(clus->Eta())>0.5) continue;
+      if((clus->EffectiveAreaCharged())<0.1*meanarea) continue; 
+      Double_t distance1=(first->Eta()-clus->Eta())*(first->Eta()-clus->Eta())+
+       (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;    
+      ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged(); 
+      count=count+1;}
+    if(count!=0) rho=ptarea/count; 
+  }        
+  return rho;
 }
 
    Double_t AliAnalysisTaskJetBackgroundSubtract::RhoRC(TClonesArray* bkgClustersRC){
@@ -586,10 +585,14 @@ Double_t AliAnalysisTaskJetBackgroundSubtract::RecalcRho(TClonesArray* bkgCluste
        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);
-        TString jetref=ostr->GetString().Data();
-        if(jetref.Contains("ANTIKT04")){ 
-        jarray = (TClonesArray*)fInJetArrayList->At(iJB);}}
+         TObjString *ostr = (TObjString*)fInJetArrayList->At(iJB);
+         TString jetref=ostr->GetString().Data();
+         if(jetref.Contains("ANTIKT04")){ 
+           jarray = (TClonesArray*)fInJetArrayList->At(iJB);
+         }
+       }
+       if(!jarray)return rho;
+
          if(jarray->GetEntries()>=2){ 
           AliAODJet *first = (AliAODJet*)(jarray->At(0)); 
           AliAODJet *second=(AliAODJet*)(jarray->At(1));