coverity fixes (R. Ma)
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 31 Jul 2013 05:01:29 +0000 (05:01 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 31 Jul 2013 05:01:29 +0000 (05:01 +0000)
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskHJetEmbed.cxx
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskHJetEmbed.h

index 524bdfa..a74933c 100644 (file)
@@ -33,6 +33,7 @@
 #include "AliGenPythiaEventHeader.h"
 #include "AliLog.h"
 #include "AliRhoParameter.h"
+#include "AliNamedString.h"
 
 #include "AliEmcalJet.h"
 
@@ -54,12 +55,13 @@ AliAnalysisTaskHJetEmbed::AliAnalysisTaskHJetEmbed() :
   fMCParticleArrName(""), fMCParticleArray(0x0), 
   fTrackArrName(""), fTrackArray(0x0), fTriggerTrkIndex(-1), 
   fMinTrkPt(0.15), fMaxTrkPt(1e4), fMinTrkEta(-0.9), fMaxTrkEta(0.9), fMinTrkPhi(0), fMaxTrkPhi(2*pi), 
-  fTTtype(0), fMinTTPt(19), fMaxTTPt(25), 
+  fTTtype(0), 
   fRadius(0.4), fJetArrName(""), fPLJetArrName(""), fDLJetArrName(""),
   fJetArray(0x0), fPLJetArray(0x0), fDLJetArray(0x0),
   fRhoName(""), fRho(0x0), fRhoValue(0), 
-  fPtHardBinParam(0), fPtHardBin(-1), fRandom(0x0),
+  fPtHardBinName(0x0), fPtHardBin(-1), fRandom(0x0),
   fRunQA(kTRUE), fRunHJet(kTRUE), fRunMatch(kTRUE),
+  fRunPL(kFALSE), fRunDL(kTRUE),
   fOutputList(0x0), fhEventStat(0x0), fhPtHardBins(0x0)
 {
   // Constructor
@@ -92,6 +94,13 @@ AliAnalysisTaskHJetEmbed::AliAnalysisTaskHJetEmbed() :
       fhJetPhiGeoMatch[i]        = 0x0;
       fhJetPhiEnMatch[i]         = 0x0;
     }
+  for(Int_t i=0; i<kNTT; i++)
+    {
+      if(i==0)
+       { fMinTTPt[i] = 19; fMaxTTPt[i] = 25; }
+      else
+       { fMinTTPt[i] = -1; fMaxTTPt[i] = -1; }
+    }
 }
 
 //________________________________________________________________________
@@ -102,12 +111,13 @@ AliAnalysisTaskHJetEmbed::AliAnalysisTaskHJetEmbed(const char *name) :
   fMCParticleArrName(""), fMCParticleArray(0x0), 
   fTrackArrName(""), fTrackArray(0x0), fTriggerTrkIndex(-1), 
   fMinTrkPt(0.15), fMaxTrkPt(1e4), fMinTrkEta(-0.9), fMaxTrkEta(0.9), fMinTrkPhi(0), fMaxTrkPhi(2*pi), 
-  fTTtype(0), fMinTTPt(19), fMaxTTPt(25), 
+  fTTtype(0), 
   fRadius(0.4), fJetArrName(""), fPLJetArrName(""), fDLJetArrName(""),
   fJetArray(0x0), fPLJetArray(0x0), fDLJetArray(0x0),
   fRhoName(""), fRho(0x0), fRhoValue(0), 
-  fPtHardBinParam(0), fPtHardBin(-1), fRandom(0x0),
+  fPtHardBinName(0x0), fPtHardBin(-1), fRandom(0x0),
   fRunQA(kTRUE), fRunHJet(kTRUE), fRunMatch(kTRUE),
+  fRunPL(kFALSE), fRunDL(kTRUE),
   fOutputList(0x0), fhEventStat(0x0), fhPtHardBins(0x0)
 {
   // Constructor
@@ -137,6 +147,13 @@ AliAnalysisTaskHJetEmbed::AliAnalysisTaskHJetEmbed(const char *name) :
       fhJetPhiGeoMatch[i]        = 0x0;
       fhJetPhiEnMatch[i]         = 0x0;
     }
+  for(Int_t i=0; i<kNTT; i++)
+    {
+      if(i==0)
+       { fMinTTPt[i] = 19; fMaxTTPt[i] = 25; }
+      else
+       { fMinTTPt[i] = -1; fMaxTTPt[i] = -1; }
+    }
 }
 //________________________________________________________________________
 AliAnalysisTaskHJetEmbed::~AliAnalysisTaskHJetEmbed()
@@ -181,10 +198,10 @@ void AliAnalysisTaskHJetEmbed::UserCreateOutputObjects()
   const Double_t lowBinMthPt[dimMthPt] = {lowJetPtBin, lowJetPtBin, -0.95, 0,  0,  0};
   const Double_t hiBinMthPt[dimMthPt]  = {upJetPtBin,  upJetPtBin,  1.05,  1,  30, 11};
 
-  const Int_t dimMthPhi = 5;
-  const Int_t nBinsMthPhi[dimMthPhi]     = {nJetPtBins,  181,                       20, 30, 11};
-  const Double_t lowBinMthPhi[dimMthPhi] = {lowJetPtBin, -pi/2-TMath::DegToRad()/2, 0,  0,  0};
-  const Double_t hiBinMthPhi[dimMthPhi]  = {upJetPtBin,  pi/2+TMath::DegToRad()/2,  1,  30, 11};
+  const Int_t dimMthPhi = 7;
+  const Int_t nBinsMthPhi[dimMthPhi]     = {nTrkPtBins,  nJetPtBins/2,  70,      70,      10, 1,  11};
+  const Double_t lowBinMthPhi[dimMthPhi] = {lowTrkPtBin, lowJetPtBin,   pi-4.95, pi-4.95, 0,  0,  0};
+  const Double_t hiBinMthPhi[dimMthPhi]  = {upTrkPtBin,  upJetPtBin,    pi+2.05, pi+2.05, 1,  10, 11};
 
   OpenFile(1);
   fOutputList = new TList();
@@ -258,10 +275,10 @@ void AliAnalysisTaskHJetEmbed::UserCreateOutputObjects()
          fhJetPtEnMatch[i] = new THnSparseF(Form("%s_fhJetPtEnMatch",triggerName[i]),Form("Embed: generated p_{T,jet} vs reconstructed p_{T,jet} vs jet p_{T} difference vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c);(p_{T,jet}^{rec}-p_{T,jet}^{gen})/p_{T,jet}^{gen};dR;centrality;pT hard bin",fRadius),dimMthPt,nBinsMthPt,lowBinMthPt,hiBinMthPt);
          fOutputList->Add(fhJetPtEnMatch[i]);
 
-         fhJetPhiGeoMatch[i] = new THnSparseF(Form("%s_fhJetPhiGeoMatch",triggerName[i]),Form("Embed: generated p_{T,jet} vs #Delta#varphi vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,jet}^{gen} (GeV/c);#Delta#varphi;dR;centrality;pT hard bin",fRadius),dimMthPhi,nBinsMthPhi,lowBinMthPhi,hiBinMthPhi);
+         fhJetPhiGeoMatch[i] = new THnSparseF(Form("%s_fhJetPhiGeoMatch",triggerName[i]),Form("Embed: p_{T,TT} vs p_{T,jet}^{det} vs #Delta#varphi_{TT} vs #Delta#varphi_{jet} vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,TT} (GeV/c);p_{T,jet}^{det} (GeV/c);#Delta#varphi_{TT};#Delta#varphi_{jet};dR;centrality;pThard bin",fRadius),dimMthPhi,nBinsMthPhi,lowBinMthPhi,hiBinMthPhi);
          fOutputList->Add(fhJetPhiGeoMatch[i]);
 
-         fhJetPhiEnMatch[i] = new THnSparseF(Form("%s_fhJetPhiEnMatch",triggerName[i]),Form("Embed: generated p_{T,jet} vs #Delta#varphi vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,jet}^{gen} (GeV/c);#Delta#varphi;dR;centrality;pT hard bin",fRadius),dimMthPhi,nBinsMthPhi,lowBinMthPhi,hiBinMthPhi);
+         fhJetPhiEnMatch[i] = new THnSparseF(Form("%s_fhJetPhiEnMatch",triggerName[i]),Form("Embed: p_{T,TT} vs p_{T,jet}^{det} vs #Delta#varphi_{TT} vs #Delta#varphi_{jet} vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,TT} (GeV/c);p_{T,jet}^{det} (GeV/c);#Delta#varphi_{TT};#Delta#varphi_{jet};dR;centrality;pThard bin",fRadius),dimMthPhi,nBinsMthPhi,lowBinMthPhi,hiBinMthPhi);
          fOutputList->Add(fhJetPhiEnMatch[i]);
        }
     }
@@ -289,31 +306,32 @@ void AliAnalysisTaskHJetEmbed::UserExec(Option_t *)
 {  
   // Main loop, called for each event.
 
+  AliDebug(5,"Entering UserExec");
   fTriggerType = -1;
-
   fEvent = InputEvent();
   if (!fEvent) 
     {
       AliError("Input event not available");
       return;
     }
-
-  if(!fPtHardBinParam)
+  AliDebug(5,"Got the input event");
+  if(!fPtHardBinName)
     {
       // Get embedded pt hard bin number
-      fPtHardBinParam = static_cast<TParameter<int>*>(fEvent->FindListObject("PYTHIAPtHardBin"));
-      if(!fPtHardBinParam)
-       {
-         AliError("The object for pt hard bin information is not available!");
-         return;
-       }
-    }
-  fPtHardBin = fPtHardBinParam->GetVal();
-  AliDebug(2,Form("Embed pt hard bin: %d\n",fPtHardBin));
-  if(fPtHardBin<0) return;
-
+      fPtHardBinName = static_cast<AliNamedString*>(fEvent->FindListObject("AODEmbeddingFile"));
+      if(!fPtHardBinName)
+       {
+         AliError("The object for pt hard bin information is not available!");
+         return;
+       }
+    }
+  TString fileName = fPtHardBinName->GetString();
+  fileName.Remove(0,50);
+  fileName.Remove(fileName.Index("/"));
+  fPtHardBin = fileName.Atoi();
   fhEventStat->Fill(0.5);
   UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+
   if(fPeriod.Contains("lhc11h",TString::kIgnoreCase))
     {
       if (trigger & AliVEvent::kAnyINT)      { fTriggerType=0; }
@@ -382,17 +400,8 @@ void AliAnalysisTaskHJetEmbed::UserExec(Option_t *)
        }
     }
 
-  const Int_t Ntracks = fTrackArray->GetEntries();
-  Int_t nLabel = 0;
-  for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) 
-    {
-      AliVParticle *t = static_cast<AliVParticle*>(fTrackArray->At(iTracks));
-      if (!t) continue;
-      if(t->GetLabel()!=0) nLabel ++;
-    }
-
   // Get MC particle array
-  if (!fMCParticleArray) 
+  if (fRunPL && !fMCParticleArray) 
     {
       fMCParticleArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fMCParticleArrName));
       if (!fMCParticleArray) 
@@ -438,7 +447,7 @@ void AliAnalysisTaskHJetEmbed::UserExec(Option_t *)
        }
     }
   // Get particle-level jet array
-  if (!fPLJetArray && !fPLJetArrName.IsNull())
+  if (fRunPL && !fPLJetArray && !fPLJetArrName.IsNull())
     {
       fPLJetArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fPLJetArrName));
       if (!fPLJetArray)
@@ -454,7 +463,7 @@ void AliAnalysisTaskHJetEmbed::UserExec(Option_t *)
        }
     }
  // Get detector-level jet array
-  if (!fDLJetArray && !fDLJetArrName.IsNull())
+  if (fRunDL && !fDLJetArray && !fDLJetArrName.IsNull())
     {
       fDLJetArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fDLJetArrName));
       if (!fDLJetArray)
@@ -475,8 +484,11 @@ void AliAnalysisTaskHJetEmbed::UserExec(Option_t *)
   fhPtHardBins->Fill(fPtHardBin);
 
   if(fRunQA) RunQA();
-  if(fRunHJet) RunHJet();
-  if(fRunMatch) RunMatch();
+  if(fRunHJet) 
+    {
+      for(Int_t i=0; i<kNTT; i++)
+       RunHJet(fMinTTPt[i],fMaxTTPt[i]);
+    }
 
   PostData(1, fOutputList);
   return;
@@ -484,50 +496,55 @@ void AliAnalysisTaskHJetEmbed::UserExec(Option_t *)
 
 
 //________________________________________________________________________
-void AliAnalysisTaskHJetEmbed::RunHJet()
+void AliAnalysisTaskHJetEmbed::RunHJet(const Double_t minPt, const Double_t maxPt)
 {
-  // Find trigger track on particle level
-  const Int_t nParticles = fMCParticleArray->GetEntries();
-  Double_t maxPLPt = -1;
-  Int_t indexPL = -1;
   TArrayI arr;
-  arr.Set(nParticles);
   Int_t counter = 0;
-  for(Int_t iPart=0; iPart<nParticles; iPart++)
+  Int_t indexPL = -1;
+
+  if(fRunPL)
     {
-      AliVParticle *t = static_cast<AliVParticle*>(fMCParticleArray->At(iPart));
-      if(!t || t->Charge()==0) continue;
-      if(!AcceptTrack(t)) continue;
-      Double_t pt = t->Pt();
-      if(fTTtype==0) // single inclusive triggers
+      // Find trigger track on particle level
+      const Int_t nParticles = fMCParticleArray->GetEntries();
+      Double_t maxPLPt = -1;
+      arr.Set(nParticles);
+      counter = 0;
+      for(Int_t iPart=0; iPart<nParticles; iPart++)
        {
-         if (pt<fMaxTTPt && pt>=fMinTTPt)
+         AliVParticle *t = static_cast<AliVParticle*>(fMCParticleArray->At(iPart));
+         //if(!t || t->Charge()==0) continue;
+         //if(!AcceptTrack(t)) continue;
+         Double_t pt = t->Pt();
+         if(fTTtype==0) // single inclusive triggers
+           {
+             if (pt<maxPt && pt>=minPt)
+               {
+                 arr.AddAt(iPart,counter);
+                 counter++;
+               }
+           }
+         else if(fTTtype==1) // leading triggers
            {
-             arr.AddAt(iPart,counter);
-             counter++;
+             if(maxPLPt<pt)
+               {
+                 maxPLPt = pt;
+                 indexPL = iPart;
+               }
            }
        }
-      else if(fTTtype==1) // leading triggers
+      arr.Set(counter);
+      if(fTTtype==0)
        {
-         if(maxPLPt<pt)
+         if(counter==0) indexPL = -1;
+         else if(counter==1) indexPL = arr.At(0);
+         else
            {
-             maxPLPt = pt;
-             indexPL = iPart;
+             Double_t pro = fRandom->Uniform() * counter;
+             indexPL = arr.At(TMath::FloorNint(pro));
            }
        }
+      arr.Reset();
     }
-  arr.Set(counter);
-  if(fTTtype==0)
-    {
-      if(counter==0) indexPL = -1;
-      else if(counter==1) indexPL = arr.At(0);
-      else
-       {
-         Double_t pro = fRandom->Uniform() * counter;
-         indexPL = arr.At(TMath::FloorNint(pro));
-       }
-    }
-  arr.Reset();
 
 
   // Find trigger track on detector level and after embedding
@@ -541,13 +558,13 @@ void AliAnalysisTaskHJetEmbed::RunHJet()
       AliVParticle *t = static_cast<AliVParticle*>(fTrackArray->At(iTracks));
       if(!t || t->Charge()==0) continue;
       if(!AcceptTrack(t)) continue;
-
       if(t->GetLabel()!=0) 
        {
+         //cout<<iTracks<<"  "<<t->Pt()<<" "<<t->GetLabel()<<endl;
          Double_t pt = t->Pt();
          if(fTTtype==0) 
            {
-             if (pt<fMaxTTPt && pt>=fMinTTPt)
+             if (pt<maxPt && pt>=minPt)
                {
                  arr.AddAt(iTracks,counter);
                  counter++;
@@ -576,12 +593,14 @@ void AliAnalysisTaskHJetEmbed::RunHJet()
     }
   arr.Reset();
 
-  AliDebug(5,Form("TT indices: PL=%d, DL=%d\n",indexPL,indexDL));
+  AliDebug(2,Form("TT indices: PL=%d, DL=%d\n",indexPL,indexDL));
 
   // Run h+jet
-  FillHJetCor(fMCParticleArray, indexPL, fPLJetArray, fhPLTT[fTriggerType], fhPLHJet[fTriggerType], kFALSE);
-  FillHJetCor(fTrackArray,      indexDL, fDLJetArray, fhDLTT[fTriggerType], fhDLHJet[fTriggerType], kFALSE);
+  if(fRunPL)  FillHJetCor(fMCParticleArray, indexPL, fPLJetArray, fhPLTT[fTriggerType], fhPLHJet[fTriggerType], kFALSE);
+  if(fRunDL)  FillHJetCor(fTrackArray,      indexDL, fDLJetArray, fhDLTT[fTriggerType], fhDLHJet[fTriggerType], kFALSE);
   FillHJetCor(fTrackArray,      indexDL, fJetArray,   fhTTPt[fTriggerType], fhHJet[fTriggerType],   kTRUE);  
+
+  if(fRunMatch) RunMatch(fTrackArray, indexDL);
 }
 
 //________________________________________________________________________
@@ -593,6 +612,7 @@ void  AliAnalysisTaskHJetEmbed::FillHJetCor(const TClonesArray *tracks, const In
   Double_t triggerPt = tt->Pt();
   Double_t fill1[] = {triggerPt, fCentrality, fPtHardBin };
   hTT->Fill(fill1);
+  AliDebug(2,Form("Found a trigger with pt = %2.2f",triggerPt));
 
   Double_t triggerPhi = tt->Phi();
   if(triggerPhi<0) triggerPhi += 2*pi;
@@ -600,46 +620,50 @@ void  AliAnalysisTaskHJetEmbed::FillHJetCor(const TClonesArray *tracks, const In
   for(Int_t ij=0; ij<nJets; ij++)
     {
       AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(jetArray->At(ij));
-      if(jet==0 || !IsGoodJet(jet)) continue; // eta cut  // jet == 0 should not occur, but Coverity wants us to check...
+      if(!jet) continue;
+      if(!IsGoodJet(jet)) continue; // eta cut
       Double_t jetPhi = jet->Phi();
       Double_t jetPt  = jet->Pt();
       Double_t jetArea = jet->Area();
       Double_t dPhi = CalculateDPhi(triggerPhi,jetPhi);
       Double_t fill[] = {triggerPt,jetPt-jetArea*fRhoValue,dPhi,jetArea,fCentrality,fPtHardBin};
       if(!isBkg) fill[1] = jetPt; 
+      AliDebug(10,"Fill the histograms");
       hn->Fill(fill);
     }
 }
 
 //________________________________________________________________________
-void AliAnalysisTaskHJetEmbed::RunMatch()
+void AliAnalysisTaskHJetEmbed::RunMatch(const TClonesArray *tracks, const Int_t leadingIndex)
 {
-  if(!fPLJetArray || !fJetArray)
+  if(leadingIndex<0) return;
+
+  if(!fDLJetArray || !fJetArray)
     {
       AliWarning("Jet array is not available.");
       return;
     }
-  Double_t dR = 999;
-  Int_t nJets = fJetArray->GetEntries();
+
+  AliVParticle *tt = (AliVParticle*) tracks->At(leadingIndex);
+  Double_t dR = 999, fraction = -1;
+  Int_t nJets = fDLJetArray->GetEntries();
   for(Int_t ij=0; ij<nJets; ij++)
     {
-      AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJetArray->At(ij));
+      AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fDLJetArray->At(ij));
+      if(!jet) continue;
       if(!IsGoodJet(jet)) continue; // eta cut
       Double_t jetPt = jet->Pt();
+      if(jetPt<10) continue;
       
-      // Geometrial matching
-      Int_t mthJetIndexGeo = FindGeoMatchedJet(jet,fPLJetArray,dR);
-      if(mthJetIndexGeo>-1)
+      // energy matching
+      Int_t mthJetIndexEn = FindEnergyMatchedJet(jet,fJetArray,dR,fraction);
+      if(mthJetIndexEn>-1 && fraction>0.5)
        {
-         AliEmcalJet* jetMthGeo = dynamic_cast<AliEmcalJet*>(fPLJetArray->At(mthJetIndexGeo));
-         Int_t dataJetIndex = FindGeoMatchedJet(jetMthGeo,fJetArray,dR);
-         if(dataJetIndex==ij) // one-to-one match
+         AliEmcalJet* jetMthEn = dynamic_cast<AliEmcalJet*>(fJetArray->At(mthJetIndexEn));
+         if(jetMthEn)
            {
-             Double_t jetPtMthGeo = jetMthGeo->Pt();
-             Double_t fill1[] = {jetPtMthGeo,jetPt,(jetPtMthGeo-jetPt)/jetPtMthGeo,dR,fCentrality, fPtHardBin};
-             fhJetPtGeoMatch[fTriggerType]->Fill(fill1);
-             Double_t fill2[] = {jetPtMthGeo,jetMthGeo->Phi()-jet->Phi(),dR,fCentrality, fPtHardBin};
-             fhJetPhiGeoMatch[fTriggerType]->Fill(fill2);
+             Double_t fill[] = {tt->Pt(),jetPt,CalculateDPhi(tt->Phi(),jet->Phi()),CalculateDPhi(jetMthEn->Phi(),jet->Phi()),dR,fCentrality,fPtHardBin};
+             fhJetPhiEnMatch[fTriggerType]->Fill(fill);
            }
        }
     }
@@ -658,6 +682,7 @@ Int_t AliAnalysisTaskHJetEmbed::FindGeoMatchedJet(const AliEmcalJet* jet, const
   for(Int_t ij=0; ij<nJets; ij++)
     {
       AliEmcalJet* jetTmp = dynamic_cast<AliEmcalJet*>(jetArray->At(ij)); 
+      if(!jetTmp) continue;
       if(TMath::Abs(jetTmp->Eta())>1) continue; // Generous eta cut
       Double_t dPhi = GetDPhi(jet->Phi(),jetTmp->Phi());
       Double_t dEta = jet->Eta()-jetTmp->Eta();
@@ -673,46 +698,51 @@ Int_t AliAnalysisTaskHJetEmbed::FindGeoMatchedJet(const AliEmcalJet* jet, const
 }
 
 //________________________________________________________________________
-Int_t AliAnalysisTaskHJetEmbed::FindEnergyMatchedJet(const AliEmcalJet* jet, const TClonesArray *jetArray, Double_t &dR)
+Int_t AliAnalysisTaskHJetEmbed::FindEnergyMatchedJet(const AliEmcalJet* jet, const TClonesArray *jetArray, Double_t &dR, Double_t &fraction)
 {
   dR = 999;
-  if(!jetArray) return -1;
+  fraction=-1;
+  if(!jetArray || !jet) return -1;
   
   Int_t index = -1;
   Int_t nJets = jetArray->GetEntries();
-  Double_t fMin = 0;
+  Double_t maxFrac = 0;
   Int_t nJetC = (Int_t)jet->GetNumberOfConstituents();
+  Double_t jetPt = jet->Pt();
   for(Int_t ij=0; ij<nJets; ij++)
     {
       AliEmcalJet* jetTmp = dynamic_cast<AliEmcalJet*>(jetArray->At(ij)); 
+      if(!jetTmp) continue;
       if(TMath::Abs(jetTmp->Eta())>1) continue; // Generous eta cut
-      Double_t jetPt = jet->Pt();
+      if(GetJetDistance(jet,jetTmp)>1) continue;
+
       Int_t nc = (Int_t)jetTmp->GetNumberOfConstituents();
       Double_t sumPt = 0;
       for(Int_t ic=0; ic<nc; ic++)
        {
-         AliVParticle *part = jetTmp->TrackAt(ic, fMCParticleArray);
          for(Int_t ijc=0; ijc<nJetC; ijc++)
            {
-             AliVParticle *track = jet->TrackAt(ijc, fTrackArray);
-             if(track->GetLabel()==part->GetLabel())
-               sumPt += part->Pt();
+             if(jetTmp->TrackAt(ic)==jet->TrackAt(ijc))
+               {
+                 AliVParticle *part = (AliVParticle*)jet->TrackAt(ijc,fTrackArray);
+                 sumPt += part->Pt();
+               }
            }
        }
       Double_t frac = sumPt/jetPt;
-      if(frac>fMin)
+      if(frac>maxFrac)
        {
-         fMin = frac;
+         maxFrac = frac;
          index = ij;
        }
     }
+  fraction = maxFrac;
 
   if(index>0)
     {
       AliEmcalJet* jetTmp = dynamic_cast<AliEmcalJet*>(jetArray->At(index)); 
-      Double_t dPhi = GetDPhi(jet->Phi(),jetTmp->Phi());
-      Double_t dEta = jet->Eta()-jetTmp->Eta();
-      dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
+      if(jetTmp)
+       dR = GetJetDistance(jet,jetTmp);
     }
   return index;
 }
@@ -738,6 +768,14 @@ Double_t AliAnalysisTaskHJetEmbed::GetDPhi(const Double_t phi1, const Double_t p
 }
 
 //________________________________________________________________________
+Double_t AliAnalysisTaskHJetEmbed::GetJetDistance(const AliEmcalJet *jet1, const AliEmcalJet* jet2)
+{
+  Double_t dPhi = GetDPhi(jet1->Phi(),jet2->Phi());
+  Double_t dEta = jet1->Eta()-jet2->Eta();
+  return TMath::Sqrt(dPhi*dPhi+dEta*dEta);
+}
+
+//________________________________________________________________________
 void AliAnalysisTaskHJetEmbed::RunQA()
 {
   if(!fPLJetArray)
@@ -750,6 +788,7 @@ void AliAnalysisTaskHJetEmbed::RunQA()
       for(Int_t ij=0; ij<nPLJets; ij++)
        {
          AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fPLJetArray->At(ij));
+         if(!jet) continue;
          if(!IsGoodJet(jet)) continue; // eta cut
          Double_t jetPt = jet->Pt();
          Double_t fill[] = {jetPt, fCentrality, fPtHardBin};
@@ -760,7 +799,7 @@ void AliAnalysisTaskHJetEmbed::RunQA()
 
   if(!fDLJetArray)
     {
-      AliWarning(Form("Particle-level jet array is not available: %s\n",fDLJetArrName.Data()));
+      AliWarning(Form("Detector-level jet array is not available: %s\n",fDLJetArrName.Data()));
     }
   else
     {
@@ -768,6 +807,7 @@ void AliAnalysisTaskHJetEmbed::RunQA()
       for(Int_t ij=0; ij<nDLJets; ij++)
        {
          AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fDLJetArray->At(ij));
+         if(!jet) continue;
          if(!IsGoodJet(jet)) continue; // eta cut
          Double_t jetPt = jet->Pt();
          Double_t fill[] = {jetPt, fCentrality, fPtHardBin};
@@ -807,8 +847,12 @@ void AliAnalysisTaskHJetEmbed::PrintConfig()
   printf("Track pt range: %2.2f < pt < %2.2f\n",fMinTrkPt, fMaxTrkPt);
   printf("Track eta range: %2.1f < eta < %2.1f\n",fMinTrkEta, fMaxTrkEta);
   printf("Track phi range: %2.0f < phi < %2.0f\n",fMinTrkPhi*TMath::RadToDeg(),fMaxTrkPhi*TMath::RadToDeg());
-  printf("TT range: %s, %2.0f < pt < %2.0f\n", TTtype[fTTtype], fMinTTPt, fMaxTTPt);
+  printf("TT type: %s\n", TTtype[fTTtype]);
+  for(Int_t i=0; i<kNTT; i++)
+    printf("TT range %d:  %2.0f < pt < %2.0f\n", i+1, fMinTTPt[i], fMaxTTPt[i]);
   printf("Run QA: %s\n",decision[fRunQA]);
+  printf("Run particle level: %s\n",decision[fRunPL]);
+  printf("Run detector level: %s\n",decision[fRunDL]);
   printf("Run h+jet: %s\n",decision[fRunHJet]);
   printf("Run matching: %s\n",decision[fRunMatch]);
   printf("=======================================\n\n");
index 3e1693b..d1055db 100644 (file)
@@ -11,6 +11,7 @@ class THnSparse;
 class TClonesArray;
 class TObject;
 class TString;
+class AliNamedString;
 class AliAODEvent;
 class AliESDEvent;
 class AliMCEvent;
@@ -42,7 +43,6 @@ class AliAnalysisTaskHJetEmbed : public AliAnalysisTaskSE {
   void SetTrkPtRange(Double_t min, Double_t max)  { fMinTrkPt=min; fMaxTrkPt=max;   }
   void SetTrkPhiRange(Double_t min, Double_t max) { fMinTrkPhi=min; fMaxTrkPhi=max; }
   void SetTrkEtaRange(Double_t min, Double_t max) { fMinTrkEta=min; fMaxTrkEta=max; }
-  void SetTTRange(Double_t min, Double_t max)     { fMinTTPt=min; fMaxTTPt=max;     }
   void SetTTtype(Int_t type)                      { fTTtype=type;                   }
   void SetRadius(Double_t rad)                    { fRadius=rad;                    }
   void SetPLJetArrName(char *s)                   { fPLJetArrName=s;                }
@@ -50,25 +50,37 @@ class AliAnalysisTaskHJetEmbed : public AliAnalysisTaskSE {
   void SetJetArrName(char *s)                     { fJetArrName=s;                  }
   void SetRhoName(char *s)                        { fRhoName=s;                     }
   void SetRunQA(Bool_t run)                       { fRunQA=run;                     }
+  void SetRunPL(Bool_t run)                       { fRunPL=run;                     }
+  void SetRunDL(Bool_t run)                       { fRunDL=run;                     }
   void SetRunHJet(Bool_t run)                     { fRunHJet=run;                   }
   void SetRunMatch(Bool_t run)                    { fRunMatch=run;                  }
+  void SetTTRanges(Double_t *min, Double_t *max)
+  {
+    for(Int_t i=0; i<kNTT; i++)
+      {
+       fMinTTPt[i] = min[i];
+       fMaxTTPt[i] = max[i];
+      }
+  }
 
 
 
 protected:
   void         RunQA();
-  void         RunHJet();
-  void         RunMatch();
+  void         RunHJet(const Double_t minPt, const Double_t maxPt);
+  void         RunMatch(const TClonesArray *tracks, const Int_t leadingIndex);
   void         FillHJetCor(const TClonesArray *tracks, const Int_t leadingIndex, const TClonesArray *jetArray, THnSparse *hTT, THnSparse *hn, Bool_t isBkg = kFALSE);
   Int_t        FindGeoMatchedJet(const AliEmcalJet* jet, const TClonesArray *jetArray, Double_t &dR);
-  Int_t        FindEnergyMatchedJet(const AliEmcalJet* jet, const TClonesArray *jetArray, Double_t &dR);
+  Int_t        FindEnergyMatchedJet(const AliEmcalJet* jet, const TClonesArray *jetArray, Double_t &dR, Double_t &fraction);
   Bool_t       AcceptTrack(const AliVParticle *track);
   Bool_t       IsGoodJet(const AliEmcalJet* jet);
   Double_t     GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz);
   Double_t     GetDPhi(const Double_t phi1, const Double_t phi2);
   Double_t     CalculateDPhi(const Double_t phi1, const Double_t phi2);
+  Double_t     GetJetDistance(const AliEmcalJet *jet1, const AliEmcalJet* jet2);
 
   enum { kNTrig = 3  };
+  enum { kNTT = 4 };
 
  private:
   Int_t             fVerbosity;            //! Control output
@@ -91,8 +103,8 @@ protected:
   Double_t          fMinTrkPhi;            //
   Double_t          fMaxTrkPhi;            //
   Int_t             fTTtype;               //
-  Double_t          fMinTTPt;              //
-  Double_t          fMaxTTPt;              //
+  Double_t          fMinTTPt[kNTT];        //
+  Double_t          fMaxTTPt[kNTT];        //
   Double_t          fRadius;               //  Jet radius
   TString           fJetArrName;           //  Name of the found jet array
   TString           fPLJetArrName;         //  Name of the embedded PYTHIA jet array on particle level
@@ -103,13 +115,15 @@ protected:
   TString           fRhoName;              //  Name of the rho parameter
   AliRhoParameter   *fRho;                 //! Rho parameter
   Double_t          fRhoValue;             //! Value of the rho parameter
-  TParameter<int>   *fPtHardBinParam;      //!Pt hard bin param
+  AliNamedString    *fPtHardBinName;       //!Pt hard bin param
   Int_t             fPtHardBin;            //!        
   TRandom3          *fRandom;              //!
 
   Bool_t            fRunQA;                //  Flag to run QA
   Bool_t            fRunHJet;              //  Flag to run h+jet 
   Bool_t            fRunMatch;             //  Flag to run matching
+  Bool_t            fRunPL;                //
+  Bool_t            fRunDL;                //
 
   TList             *fOutputList;                //! Output list
   TH1F              *fhEventStat;                //!
@@ -143,7 +157,9 @@ protected:
   AliAnalysisTaskHJetEmbed(const AliAnalysisTaskHJetEmbed&);            // not implemented
   AliAnalysisTaskHJetEmbed &operator=(const AliAnalysisTaskHJetEmbed&); // not implemented
 
-  ClassDef(AliAnalysisTaskHJetEmbed, 2);
+  ClassDef(AliAnalysisTaskHJetEmbed, 3);
 };
 
 #endif
+
+