bug-fix: rotation of sub-leading jet in di-jet

[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskEmcalJetTagger.cxx

diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetTagger.cxx b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetTagger.cxx
index f0fea8a51a465525caf1531b76236a8c3bff3b36..47b4efd25f479dc5128ac1c8add971d43c445d91 100644 (file)
--- a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetTagger.cxx
+++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetTagger.cxx
@@ -52,6 +52,7 @@ AliAnalysisTaskEmcalJetTagger::AliAnalysisTaskEmcalJetTagger() :
   fh2PtJet1VsLeadPtAllSel(0),
   fh2PtJet1VsLeadPtTagged(0),
   fh2PtJet1VsPtJet2(0),
+  fh2PtJet2VsRelPt(0),
   fh3PtJetDEtaDPhiConst(0),
   fh2PtJetDRConst(0),
   fh3PtJetAreaDRConst(0)
@@ -64,6 +65,7 @@ AliAnalysisTaskEmcalJetTagger::AliAnalysisTaskEmcalJetTagger() :
   fh2PtJet1VsLeadPtAllSel      = new TH2F*[fNcentBins];
   fh2PtJet1VsLeadPtTagged      = new TH2F*[fNcentBins];
   fh2PtJet1VsPtJet2            = new TH2F*[fNcentBins];
+  fh2PtJet2VsRelPt             = new TH2F*[fNcentBins];
 
   for (Int_t i = 0; i < fNcentBins; i++) {
     fh3PtJet1VsDeltaEtaDeltaPhi[i] = 0;
@@ -72,10 +74,10 @@ AliAnalysisTaskEmcalJetTagger::AliAnalysisTaskEmcalJetTagger() :
     fh2PtJet1VsLeadPtAllSel[i]     = 0;
     fh2PtJet1VsLeadPtTagged[i]     = 0;
     fh2PtJet1VsPtJet2[i]           = 0;
+    fh2PtJet2VsRelPt[i]            = 0;
   }
 
   SetMakeGeneralHistograms(kTRUE);
-  
 }
 
 //________________________________________________________________________
@@ -93,6 +95,7 @@ AliAnalysisTaskEmcalJetTagger::AliAnalysisTaskEmcalJetTagger(const char *name) :
   fh2PtJet1VsLeadPtAllSel(0),
   fh2PtJet1VsLeadPtTagged(0),
   fh2PtJet1VsPtJet2(0),
+  fh2PtJet2VsRelPt(0),
   fh3PtJetDEtaDPhiConst(0),
   fh2PtJetDRConst(0),
   fh3PtJetAreaDRConst(0)
@@ -105,6 +108,7 @@ AliAnalysisTaskEmcalJetTagger::AliAnalysisTaskEmcalJetTagger(const char *name) :
   fh2PtJet1VsLeadPtAllSel     = new TH2F*[fNcentBins];
   fh2PtJet1VsLeadPtTagged     = new TH2F*[fNcentBins];
   fh2PtJet1VsPtJet2           = new TH2F*[fNcentBins];
+  fh2PtJet2VsRelPt            = new TH2F*[fNcentBins];
 
   for (Int_t i = 0; i < fNcentBins; i++) {
     fh3PtJet1VsDeltaEtaDeltaPhi[i] = 0;
@@ -113,6 +117,7 @@ AliAnalysisTaskEmcalJetTagger::AliAnalysisTaskEmcalJetTagger(const char *name) :
     fh2PtJet1VsLeadPtAllSel[i]     = 0;
     fh2PtJet1VsLeadPtTagged[i]     = 0;
     fh2PtJet1VsPtJet2[i]           = 0;
+    fh2PtJet2VsRelPt[i]            = 0;
   }
 
   SetMakeGeneralHistograms(kTRUE);
@@ -185,6 +190,10 @@ void AliAnalysisTaskEmcalJetTagger::UserCreateOutputObjects()
     fh2PtJet1VsPtJet2[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
     fOutput->Add(fh2PtJet1VsPtJet2[i]);
 
+    histName = TString::Format("fh2PtJet2VsRelPt_%d",i);
+    histTitle = TString::Format("%s;#it{p}_{T,jet2};(#it{p}_{T,jet2}-#it{p}_{T,jet1})/#it{p}_{T,jet1}",histName.Data());
+    fh2PtJet2VsRelPt[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,241,-2.41,2.41);
+    fOutput->Add(fh2PtJet2VsRelPt[i]);
   }
 
   fh3PtJetDEtaDPhiConst = new TH3F("fh3PtJetDEtaDPhiConst","fh3PtJetDEtaDPhiConst;pT;#Delta #eta;#Delta #varphi",nBinsPt,minPt,maxPt,nBinsDEta,-1.,1.,nBinsDPhi,-1.,1.);
@@ -217,7 +226,7 @@ Bool_t AliAnalysisTaskEmcalJetTagger::Run()
 {
   // Run analysis code here, if needed. It will be executed before FillHistograms().
 
-  MatchJetsGeo(fContainerBase,fContainerTag,0,0.3,2);
+  MatchJetsGeo(fContainerBase,fContainerTag,0,0.3,3);
 
   //  if(fJetTaggingMethod==kFraction)
 
@@ -269,6 +278,7 @@ Bool_t AliAnalysisTaskEmcalJetTagger::FillHistograms()
 
     fh2PtJet1VsLeadPtTagged[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
     fh2PtJet1VsPtJet2[fCentBin]->Fill(ptJet1,ptJet2);
+    if(ptJet2>0.) fh2PtJet2VsRelPt[fCentBin]->Fill(ptJet2,(ptJet1-ptJet2)/ptJet2);
 
     Double_t dPhi = GetDeltaPhi(jet1->Phi(),jet2->Phi());
     if(dPhi>TMath::Pi())
@@ -305,7 +315,11 @@ void AliAnalysisTaskEmcalJetTagger::MatchJetsGeo(Int_t c1, Int_t c2,
   //
   // Match the full jets to the corresponding charged jets
   // Translation of AliAnalysisHelperJetTasks::GetClosestJets to AliEmcalJet objects
-  // type: 0 = use acceptance cuts of container  1 = allow 0.1 one more for c2 in eta 2 = allow 0.1 more in eta and phi for c2
+  // type: 
+  //         0 = use acceptance cuts of container  
+  //         1 = allow 0.1 one more for c2 in eta 
+  //         2 = allow 0.1 more in eta and phi for c2
+  //         3 = allow 0.1 in eta and phi for both containers
 
   if(c1<0) c1 = fContainerBase;
   if(c2<0) c2 = fContainerTag;
@@ -331,11 +345,25 @@ void AliAnalysisTaskEmcalJetTagger::MatchJetsGeo(Int_t c1, Int_t c2,
 
   static TArrayS iFlag(nJets1*nJets2);
   if(iFlag.GetSize()<(nJets1*nJets2)){
-    iFlag.Set(nJets1*nJets1+1);
+    iFlag.Set(nJets1*nJets2+1);
   }
   iFlag.Reset(0);
 
+
+  AliJetContainer *cont1 = GetJetContainer(c1);
   AliJetContainer *cont2 = GetJetContainer(c2);
+  if(type==1)
+    cont2->SetJetEtaLimits(cont2->GetJetEtaMin()-0.1,cont2->GetJetEtaMax()+0.1);
+  else if(type==2) {
+    cont2->SetJetEtaLimits(cont2->GetJetEtaMin()-0.1,cont2->GetJetEtaMax()+0.1);
+    cont2->SetJetPhiLimits(cont2->GetJetPhiMin()-0.1,cont2->GetJetPhiMax()+0.1);
+  } 
+  else if(type==3) {
+    cont1->SetJetEtaLimits(cont1->GetJetEtaMin()-0.1,cont1->GetJetEtaMax()+0.1);
+    cont1->SetJetPhiLimits(cont1->GetJetPhiMin()-0.1,cont1->GetJetPhiMax()+0.1);
+    cont2->SetJetEtaLimits(cont2->GetJetEtaMin()-0.1,cont2->GetJetEtaMax()+0.1);
+    cont2->SetJetPhiLimits(cont2->GetJetPhiMin()-0.1,cont2->GetJetPhiMax()+0.1);
+  }
 
   // find the closest distance to the full jet
   for(int i = 0;i<nJets1;i++){
@@ -346,25 +374,10 @@ void AliAnalysisTaskEmcalJetTagger::MatchJetsGeo(Int_t c1, Int_t c2,
     Float_t dist = maxDist;
     
     for(int j = 0;j <nJets2; j++){
-      AliEmcalJet *jet2 = 0x0;
-      if(type==0)
-       jet2 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(j, c2));
-      else {
-       jet2 = static_cast<AliEmcalJet*>(GetJetFromArray(j, c2));
-       if(!jet2) continue;
-       if(type>0) {
-         if(jet2->Eta()<(cont2->GetJetEtaMin()-0.1) || jet2->Eta()>(cont2->GetJetEtaMax()+0.1))
-           continue;
-         if(type==2) {
-           if(jet2->Phi()<(cont2->GetJetPhiMin()-0.1) || jet2->Phi()>(cont2->GetJetPhiMax()+0.1))
-             continue;
-         }
-       }
-      }
-      if(!jet2)
-       continue;
+      AliEmcalJet *jet2 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(j, c2));
+      if(!jet2) continue;
 
-      Double_t dR = GetDeltaR(jet1,jet2);
+      Double_t dR = jet1->DeltaR(jet2);//GetDeltaR(jet1,jet2);
       if(dR<dist && dR<maxDist){
        faMatchIndex2[i]=j;
        dist = dR;
@@ -380,31 +393,16 @@ void AliAnalysisTaskEmcalJetTagger::MatchJetsGeo(Int_t c1, Int_t c2,
 
   // other way around
   for(int j = 0;j<nJets2;j++){
-    AliEmcalJet *jet2 = 0x0;
-    if(type==0)
-      jet2 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(j, c2));
-    else {
-      jet2 = static_cast<AliEmcalJet*>(GetJetFromArray(j, c2));
-      if(!jet2) continue;;
-      if(type>0) {
-       if(jet2->Eta()<(cont2->GetJetEtaMin()-0.1) || jet2->Eta()>(cont2->GetJetEtaMax()+0.1))
-         continue;
-       if(type==2) {
-         if(jet2->Phi()<(cont2->GetJetPhiMin()-0.1) || jet2->Phi()>(cont2->GetJetPhiMax()+0.1))
-           continue;
-       }
-      }
-    }
+    AliEmcalJet *jet2 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(j, c2));
     if(!jet2)
       continue;
 
     Float_t dist = maxDist;
     for(int i = 0;i<nJets1;i++){
       AliEmcalJet *jet1 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, c1));
-      if(!jet1)
-       continue;
+      if(!jet1)        continue;
 
-      Double_t dR = GetDeltaR(jet1,jet2); 
+      Double_t dR = jet1->DeltaR(jet2);//GetDeltaR(jet1,jet2); 
       if(dR<dist && dR<maxDist){
        faMatchIndex1[j]=i;
         dist = dR;
@@ -427,7 +425,7 @@ void AliAnalysisTaskEmcalJetTagger::MatchJetsGeo(Int_t c1, Int_t c2,
       // we have a uniqe correlation
       if(iFlag[i*nJets2+j]==3){
        
-       Double_t dR = GetDeltaR(jet1,jet2); 
+       Double_t dR = jet1->DeltaR(jet2);//GetDeltaR(jet1,jet2); 
        if(iDebug>1) Printf("closest jets %d  %d  dR =  %f",j,i,dR);
        
        if(fJetTaggingType==kTag) {