X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=PWGJE%2FEMCALJetTasks%2FAliEmcalJetTask.cxx;h=dddc0c76ee228e1b56510608362e2ede998b0d6d;hb=bebb68db31262a7b6712284a98a9500d32ff64f1;hp=d75eeacb68399bc5e98f69504a7de980d9e2fc53;hpb=816bb2c86d38478fdf656d0a19b1de6c2d1b6ec5;p=u%2Fmrichter%2FAliRoot.git

diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx b/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx
index d75eeacb683..dddc0c76ee2 100644
--- a/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx
+++ b/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx
@@ -1,8 +1,7 @@
-// $Id$
 //
 // Emcal jet finder task.
 //
-// Authors: C.Loizides, S.Aiola
+// Authors: C.Loizides, S.Aiola, M. Verweij
 
 #include <vector>
 #include "AliEmcalJetTask.h"
@@ -25,6 +24,7 @@
 #include "AliVCluster.h"
 #include "AliVEvent.h"
 #include "AliVParticle.h"
+#include "AliRhoParameter.h"
 using std::cout;
 using std::endl;
 using std::cerr;
@@ -37,6 +37,7 @@ AliEmcalJetTask::AliEmcalJetTask() :
   fTracksName("Tracks"),
   fCaloName("CaloClusters"),
   fJetsName("Jets"),
+  fJetsSubName(""),
   fJetType(kNone),
   fConstSel(0),
   fMCConstSel(0),
@@ -58,16 +59,33 @@ AliEmcalJetTask::AliEmcalJetTask() :
   fMinMCLabel(0),
   fRecombScheme(fastjet::pt_scheme),
   fTrackEfficiency(1.),
+  fMCFlag(0),
+  fGeneratorIndex(-1),
   fIsInit(0),
+  fLocked(0),
   fIsPSelSet(0),
   fIsMcPart(0),
   fIsEmcPart(0),
   fLegacyMode(kFALSE),
   fCodeDebug(kFALSE),
+  fDoGenericSubtractionJetMass(kFALSE),
+  fDoGenericSubtractionGR(kFALSE),
+  fDoConstituentSubtraction(kFALSE),
+  fUseExternalBkg(kFALSE),
+  fRhoName(""),
+  fRhomName(""),
+  fRho(0),
+  fRhom(0),
+  fRMax(0.4),
+  fDRStep(0.04),
+  fPtMinGR(40.),
   fJets(0),
+  fJetsSub(0),
   fEvent(0),
   fTracks(0),
-  fClus(0)
+  fClus(0),
+  fRhoParam(0),
+  fRhomParam(0)
 {
   // Default constructor.
 }
@@ -78,6 +96,7 @@ AliEmcalJetTask::AliEmcalJetTask(const char *name) :
   fTracksName("Tracks"),
   fCaloName("CaloClusters"),
   fJetsName("Jets"),
+  fJetsSubName(""),
   fJetType(kAKT|kFullJet|kRX1Jet),
   fConstSel(0),
   fMCConstSel(0),
@@ -99,16 +118,33 @@ AliEmcalJetTask::AliEmcalJetTask(const char *name) :
   fMinMCLabel(0),
   fRecombScheme(fastjet::pt_scheme),
   fTrackEfficiency(1.),
+  fMCFlag(0),
+  fGeneratorIndex(-1),
   fIsInit(0),
+  fLocked(0),
   fIsPSelSet(0),
   fIsMcPart(0),
   fIsEmcPart(0),
   fLegacyMode(kFALSE),
   fCodeDebug(kFALSE),
+  fDoGenericSubtractionJetMass(kFALSE),
+  fDoGenericSubtractionGR(kFALSE),
+  fDoConstituentSubtraction(kFALSE),
+  fUseExternalBkg(kFALSE),
+  fRhoName(""),
+  fRhomName(""),
+  fRho(0),
+  fRhom(0),
+  fRMax(0.4),
+  fDRStep(0.04),
+  fPtMinGR(40.),
   fJets(0),
+  fJetsSub(0),
   fEvent(0),
   fTracks(0),
-  fClus(0)
+  fClus(0),
+  fRhoParam(0),
+  fRhomParam(0)
 {
   // Standard constructor.
 
@@ -128,6 +164,11 @@ void AliEmcalJetTask::UserCreateOutputObjects()
 
   fJets = new TClonesArray("AliEmcalJet");
   fJets->SetName(fJetsName);
+
+  if(!fJetsSubName.IsNull()) {
+    fJetsSub = new TClonesArray("AliEmcalJet");
+    fJetsSub->SetName(fJetsSubName);
+  }
 }
 
 //________________________________________________________________________
@@ -142,6 +183,7 @@ void AliEmcalJetTask::UserExec(Option_t *)
 
   // clear the jet array (normally a null operation)
   fJets->Delete();
+  if(fJetsSub) fJetsSub->Delete();
 
   FindJets();
 }
@@ -161,6 +203,11 @@ void AliEmcalJetTask::FindJets()
     return;
   }
 
+ if (fRhoParam)
+    fRho = fRhoParam->GetVal();
+  if (fRhomParam)
+    fRhom = fRhomParam->GetVal();
+
   TString name("kt");
   fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
   if ((fJetType & kAKT) != 0) {
@@ -228,7 +275,15 @@ void AliEmcalJetTask::FindJets()
 	else {
 	  AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));	    
 	}
-      }	    
+      }
+      if ((t->GetFlag() & fMCFlag) != fMCFlag) { 
+	AliDebug(2,Form("Skipping track %d because it does not match the MC flags", iTracks));
+	continue;
+      }
+      if (fGeneratorIndex >= 0 && t->GetGeneratorIndex() != fGeneratorIndex) { 
+	AliDebug(2,Form("Skipping track %d because it does not match the MC generator index", iTracks));
+	continue;
+      }
       if (t->Pt() < fMinJetTrackPt) 
         continue;
       Double_t eta = t->Eta();
@@ -346,6 +401,18 @@ void AliEmcalJetTask::FindJets()
   // run jet finder
   fjw.Run();
 
+  //run generic subtractor
+  if(fDoGenericSubtractionJetMass) {
+    fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
+    fjw.DoGenericSubtractionJetMass();
+  }
+
+  //run constituent subtractor
+  if(fDoConstituentSubtraction) {
+    fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
+    fjw.DoConstituentSubtraction();
+  }
+
   // get geometry
   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
   if (!geom) {
@@ -374,6 +441,48 @@ void AliEmcalJetTask::FindJets()
 
     AliEmcalJet *jet = new ((*fJets)[jetCount]) 
       AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
+    jet->SetLabel(ij);
+
+    //do generic subtraction if requested
+#ifdef FASTJET_VERSION
+    if(fDoGenericSubtractionJetMass) {
+      std::vector<fastjet::contrib::GenericSubtractorInfo> jetMassInfo = fjw.GetGenSubtractorInfoJetMass();
+      Int_t n = (Int_t)jetMassInfo.size();
+      if(n>ij && n>0) {
+        jet->SetFirstDerivative(jetMassInfo[ij].first_derivative());
+        jet->SetSecondDerivative(jetMassInfo[ij].second_derivative());
+        jet->SetFirstOrderSubtracted(jetMassInfo[ij].first_order_subtracted());
+        jet->SetSecondOrderSubtracted(jetMassInfo[ij].second_order_subtracted());
+      }
+    }
+    //here do generic subtraction for angular structure function
+    Double_t ptcorr = jets_incl[ij].perp()-fjw.GetJetArea(ij)*fRho;
+    fRMax = fRadius+0.2;
+    if(fDoGenericSubtractionGR && ptcorr>fPtMinGR) {
+      fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
+      fjw.SetRMaxAndStep(fRMax,fDRStep);
+      fjw.DoGenericSubtractionGR(ij);
+      std::vector<double> num = fjw.GetGRNumerator();
+      std::vector<double> den = fjw.GetGRDenominator();
+      std::vector<double> nums = fjw.GetGRNumeratorSub();
+      std::vector<double> dens = fjw.GetGRDenominatorSub();
+      //pass this to AliEmcalJet
+      jet->SetGRNumSize(num.size());
+      jet->SetGRDenSize(den.size());
+      jet->SetGRNumSubSize(nums.size());
+      jet->SetGRDenSubSize(dens.size());
+      Int_t nsize = (Int_t)num.size();
+      for(Int_t g = 0; g<nsize; ++g) {
+        jet->AddGRNumAt(num[g],g);
+        jet->AddGRNumSubAt(nums[g],g);
+      }
+      Int_t dsize = (Int_t)den.size();
+      for(Int_t g = 0; g<dsize; ++g) {
+        jet->AddGRDenAt(den[g],g);
+        jet->AddGRDenSubAt(dens[g],g);
+      }
+    }
+#endif
 
     // loop over constituents
     std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
@@ -531,6 +640,30 @@ void AliEmcalJetTask::FindJets()
     jetCount++;
   }
   //fJets->Sort();
+
+  //run constituent subtractor if requested
+  if(fDoConstituentSubtraction) {
+#ifdef FASTJET_VERSION
+    if(!fJetsSub) AliWarning(Form("No jet branch to write to for subtracted jets. fJetsSubName: %s",fJetsSubName.Data()));
+    else {
+      std::vector<fastjet::PseudoJet> jets_sub;
+      jets_sub = fjw.GetConstituentSubtrJets();
+      AliDebug(1,Form("%d constituent subtracted jets found", (Int_t)jets_sub.size()));
+      for (UInt_t ijet=0, jetCount=0; ijet<jets_sub.size(); ++ijet) {
+	//Only storing 4-vector and jet area of unsubtracted jet	
+	AliEmcalJet *jet_sub = new ((*fJetsSub)[jetCount]) 
+	  AliEmcalJet(jets_sub[ijet].perp(), jets_sub[ijet].eta(), jets_sub[ijet].phi(), jets_sub[ijet].m());
+	jet_sub->SetLabel(ijet);
+	fastjet::PseudoJet area(fjw.GetJetAreaVector(ijet));
+	jet_sub->SetArea(area.perp());
+	jet_sub->SetAreaEta(area.eta());
+	jet_sub->SetAreaPhi(area.phi());
+	jet_sub->SetAreaEmc(area.perp());
+	jetCount++;
+      }
+    }
+#endif
+  } //constituent subtraction
 }
 
 //________________________________________________________________________
@@ -581,6 +714,9 @@ Bool_t AliEmcalJetTask::DoInit()
     return 0;
   }
 
+  if (!(fEvent->FindListObject(fJetsSubName)) && fJetsSub)
+    fEvent->AddObject(fJetsSub);
+
   if (fTracksName == "Tracks")
     am->LoadBranch("Tracks");
   if (!fTracks && !fTracksName.IsNull()) {
@@ -611,5 +747,20 @@ Bool_t AliEmcalJetTask::DoInit()
       fIsEmcPart = 1;
   }
 
+  if (!fRhoName.IsNull() && !fRhoParam) { // get rho from the event
+    fRhoParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
+    if (!fRhoParam) {
+      AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
+      return 0;
+    }
+  }
+  if (!fRhomName.IsNull() && !fRhomParam) { // get rhom from the event
+    fRhomParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhomName));
+    if (!fRhomParam) {
+      AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhomName.Data()));
+      return 0;
+    }
+  }
+  
   return 1;
 }