pileup rejection for Fragmentation Function and Jet Shape task
authorobusch <Oliver.Busch@cern.ch>
Thu, 17 Apr 2014 07:43:32 +0000 (09:43 +0200)
committermvl <marco.van.leeuwen@cern.ch>
Thu, 17 Apr 2014 09:12:00 +0000 (11:12 +0200)
PWGJE/AliAnalysisTaskFragmentationFunction.cxx
PWGJE/AliAnalysisTaskFragmentationFunction.h
PWGJE/AliAnalysisTaskJetProperties.cxx
PWGJE/AliAnalysisTaskJetProperties.h
PWGJE/macros/AddTaskJetProperties.C

index 1ed33b5..c3a89a4 100644 (file)
@@ -78,6 +78,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction()
    ,fEvtSelectionMask(0)
    ,fEventClass(0)
    ,fMaxVertexZ(10)
+   ,fRejectPileup(kFALSE)
    ,fTrackPtCut(0)
    ,fTrackEtaMin(0)
    ,fTrackEtaMax(0)
@@ -278,6 +279,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const
   ,fEvtSelectionMask(0)
   ,fEventClass(0)
   ,fMaxVertexZ(10)
+  ,fRejectPileup(kFALSE)
   ,fTrackPtCut(0)
   ,fTrackEtaMin(0)
   ,fTrackEtaMax(0)
@@ -479,6 +481,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const
   ,fEvtSelectionMask(copy.fEvtSelectionMask)
   ,fEventClass(copy.fEventClass)
   ,fMaxVertexZ(copy.fMaxVertexZ)
+  ,fRejectPileup(copy.fRejectPileup)
   ,fTrackPtCut(copy.fTrackPtCut)
   ,fTrackEtaMin(copy.fTrackEtaMin)
   ,fTrackEtaMax(copy.fTrackEtaMax)
@@ -684,6 +687,7 @@ AliAnalysisTaskFragmentationFunction& AliAnalysisTaskFragmentationFunction::oper
     fEvtSelectionMask              = o.fEvtSelectionMask;
     fEventClass                    = o.fEventClass;
     fMaxVertexZ                    = o.fMaxVertexZ;
+    fRejectPileup                  = o.fRejectPileup;
     fTrackPtCut                    = o.fTrackPtCut;
     fTrackEtaMin                   = o.fTrackEtaMin;
     fTrackEtaMax                   = o.fTrackEtaMax;
@@ -1396,14 +1400,15 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
   
   
   // Histograms        
-  fh1EvtSelection            = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
+  fh1EvtSelection            = new TH1F("fh1EvtSelection", "Event Selection", 7, -0.5, 6.5);
   fh1EvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
   fh1EvtSelection->GetXaxis()->SetBinLabel(2,"event selection: rejected");
   fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected");
   fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected");
   fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected");
   fh1EvtSelection->GetXaxis()->SetBinLabel(6,"vertex type: rejected");
-  
+  fh1EvtSelection->GetXaxis()->SetBinLabel(7,"pileup: rejected");
   fh1VertexNContributors     = new TH1F("fh1VertexNContributors", "Vertex N contributors", 2500,-.5, 2499.5);
   fh1VertexZ                 = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
   fh1EvtMult                = new TH1F("fh1EvtMult","Event multiplicity, track pT cut > 150 MeV/c, |#eta| < 0.9",120,0.,12000.);
@@ -2183,6 +2188,13 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *)
     return;
   }
 
+  if(fRejectPileup && AliAnalysisHelperJetTasks::IsPileUp()){
+    if (fDebug > 1) Printf("%s:%d SPD pileup: event REJECTED...",(char*)__FILE__,__LINE__);
+    fh1EvtSelection->Fill(6.);
+    PostData(1, fCommonHistList);
+    return;
+  }
+
   if (fDebug > 1) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__); 
   fh1EvtSelection->Fill(0.);
   fh1EvtCent->Fill(centPercent);
@@ -2703,6 +2715,7 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *)
          TList* perpjettracklistGen1 = new TList();
          TList* perpjettracklistGen2 = new TList();
 
+         Double_t sumPtGenPerp  = 0.;
          Double_t sumPtGenPerp1 = 0.;
          Double_t sumPtGenPerp2 = 0.;
          GetTracksTiltedwrpJetAxis(TMath::Pi()/2.,fTracksAODMCCharged, perpjettracklistGen1, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerp1); 
@@ -2710,31 +2723,36 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *)
 
          perpjettracklistGen->AddAll(perpjettracklistGen1);
          perpjettracklistGen->AddAll(perpjettracklistGen2);
+         sumPtGenPerp = 0.5*(sumPtGenPerp1+sumPtGenPerp2);
 
          TList* perpjettracklistGenSecNS  = new TList();
          TList* perpjettracklistGenSecNS1 = new TList();
          TList* perpjettracklistGenSecNS2 = new TList();
 
-          Double_t sumPtGenPerpNS1 = 0;
-          Double_t sumPtGenPerpNS2 = 0;
+          Double_t sumPtGenPerpNS;
+          Double_t sumPtGenPerpNS1;
+          Double_t sumPtGenPerpNS2;
           GetTracksTiltedwrpJetAxis(TMath::Pi()/2.,fTracksAODMCChargedSecNS, perpjettracklistGenSecNS1, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerpNS1); 
           GetTracksTiltedwrpJetAxis(-1*TMath::Pi()/2.,fTracksAODMCChargedSecNS, perpjettracklistGenSecNS2, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerpNS2); 
 
          perpjettracklistGenSecNS->AddAll(perpjettracklistGenSecNS1);
          perpjettracklistGenSecNS->AddAll(perpjettracklistGenSecNS2);
+         sumPtGenPerpNS = 0.5*(sumPtGenPerpNS1+sumPtGenPerpNS2);
 
 
          TList* perpjettracklistGenSecS  = new TList();
          TList* perpjettracklistGenSecS1 = new TList();
          TList* perpjettracklistGenSecS2 = new TList();
 
-          Double_t sumPtGenPerpS1 = 0;
-          Double_t sumPtGenPerpS2 = 0;
+          Double_t sumPtGenPerpS;
+          Double_t sumPtGenPerpS1;
+          Double_t sumPtGenPerpS2;
           GetTracksTiltedwrpJetAxis(TMath::Pi()/2.,fTracksAODMCChargedSecS, perpjettracklistGenSecS1, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerpS1); 
           GetTracksTiltedwrpJetAxis(-1*TMath::Pi()/2.,fTracksAODMCChargedSecS, perpjettracklistGenSecS2, jet, TMath::Abs(GetFFRadius()) , sumPtGenPerpS2); 
 
          perpjettracklistGenSecS->AddAll(perpjettracklistGenSecS1);
          perpjettracklistGenSecS->AddAll(perpjettracklistGenSecS2);
+         sumPtGenPerpS = 0.5*(sumPtGenPerpS1+sumPtGenPerpS2);
 
 
           if(perpjettracklistGen->GetSize() != perpjettracklistGen1->GetSize() + perpjettracklistGen2->GetSize()){
@@ -3268,7 +3286,7 @@ void AliAnalysisTaskFragmentationFunction::SetProperties(TH1* h,const char* x, c
 
 // ________________________________________________________________________________________________________________________________________________________
 void AliAnalysisTaskFragmentationFunction::GetJetTracksPointing(TList* inputlist, TList* outputlist, const AliAODJet* jet, 
-                                                               const Double_t& radius, Double_t& sumPt, const Double_t& minPtL, const Double_t& maxPt, Bool_t& isBadPt)
+                                                                  const Double_t radius, Double_t& sumPt, const Double_t minPtL, const Double_t maxPt, Bool_t& isBadPt)
 {
   // fill list of tracks in cone around jet axis  
 
@@ -3309,8 +3327,7 @@ void AliAnalysisTaskFragmentationFunction::GetJetTracksPointing(TList* inputlist
 }
 
 // _________________________________________________________________________________________________________________________________________________________________
-void AliAnalysisTaskFragmentationFunction::GetJetTracksTrackrefs(TList* list, const AliAODJet* jet, const Double_t& minPtL, 
-                                                                const Double_t& maxPt, Bool_t& isBadPt)
+void AliAnalysisTaskFragmentationFunction::GetJetTracksTrackrefs(TList* list, const AliAODJet* jet, const Double_t minPtL, const Double_t maxPt, Bool_t& isBadPt)
 {
   // list of jet tracks from trackrefs
   
@@ -3991,7 +4008,7 @@ void AliAnalysisTaskFragmentationFunction::GetClusterTracksMedian(TList* outputl
   // get median cluster
 
   AliAODJet* medianCluster = 0;
-  //Double_t   medianDensity = 0;
+  Double_t   medianDensity = 0;
 
   if(TMath::Odd(nBckgClusters)){
     
@@ -4001,7 +4018,7 @@ void AliAnalysisTaskFragmentationFunction::GetClusterTracksMedian(TList* outputl
     Double_t clusterPt = medianCluster->Pt();
     Double_t area      = medianCluster->EffectiveAreaCharged();
     
-    //if(area>0) medianDensity = clusterPt/area;
+    if(area>0) medianDensity = clusterPt/area;
   }
   else{
 
@@ -4021,7 +4038,7 @@ void AliAnalysisTaskFragmentationFunction::GetClusterTracksMedian(TList* outputl
     Double_t area2      = medianCluster2->EffectiveAreaCharged();
     if(area2>0) density2 = clusterPt2/area2;
     
-    //medianDensity = 0.5*(density1+density2);
+    medianDensity = 0.5*(density1+density2);
     
     medianCluster = ( (fRandom->Rndm()>0.5) ? medianCluster1 : medianCluster2 );  // select one randomly to avoid adding areas
   }
@@ -4614,7 +4631,7 @@ void AliAnalysisTaskFragmentationFunction::FillBckgHistos(Int_t type, TList* inp
 }
 
 //_____________________________________________________________________________________
-Double_t AliAnalysisTaskFragmentationFunction::GetMCStrangenessFactor(const Double_t& pt)
+Double_t AliAnalysisTaskFragmentationFunction::GetMCStrangenessFactor(const Double_t pt)
 {
   // factor strangeness data/MC as function of pt from UE analysis (Sara Vallero)
 
index b538e39..7bc6a80 100644 (file)
@@ -206,6 +206,7 @@ class AliAnalysisTaskFragmentationFunction : public AliAnalysisTaskSE {
   virtual void   SetEventSelectionMask(UInt_t i){fEvtSelectionMask = i;}
   virtual void   SetEventClass(Int_t i){fEventClass = i;}
   virtual void   SetMaxVertexZ(Float_t z){fMaxVertexZ = z;}
+  virtual void   RejectPileupEvents(Bool_t b){fRejectPileup = b;}
   virtual void   SetJetCuts(Float_t jetPt = 5., Float_t jetEtaMin = -0.5, Float_t jetEtaMax = 0.5, 
                            Float_t jetPhiMin = 0., Float_t jetPhiMax = 2*TMath::Pi())
   {fJetPtCut = jetPt; fJetEtaMin = jetEtaMin; fJetEtaMax = jetEtaMax; 
@@ -260,9 +261,8 @@ class AliAnalysisTaskFragmentationFunction : public AliAnalysisTaskSE {
   Float_t  GetFFMaxTrackPt() const { return fFFMaxTrackPt; }
   Float_t  GetFFMinNTracks() const { return fFFMinnTracks; }
   Float_t  GetFFBckgRadius() const { return fFFBckgRadius; }
-  void    GetJetTracksTrackrefs(TList* l, const AliAODJet* j, const Double_t& minPtL, const Double_t& maxPt, Bool_t& isBadPt);
-  void    GetJetTracksPointing(TList* in, TList* out, const AliAODJet* j, const Double_t& r, Double_t& sumPt, 
-                               const Double_t& minPtL, const Double_t& maxPt, Bool_t& isBadPt);  
+  void    GetJetTracksTrackrefs(TList* l, const AliAODJet* j, const Double_t minPtL, const Double_t maxPt, Bool_t& isBadPt);
+  void    GetJetTracksPointing(TList* in, TList* out, const AliAODJet* j, const Double_t r, Double_t& sumPt, const Double_t minPtL, const Double_t maxPt, Bool_t& isBadPt);  
   void     GetTracksOutOfNJets(Int_t nCases, TList* in, TList* out, TList* jets, Double_t& pt);
   void     GetTracksOutOfNJetsStat(Int_t nCases, TList* in, TList* out, TList* jets, Double_t& pt, Double_t &normFactor);
   void     GetTracksTiltedwrpJetAxis(Float_t alpha, TList* inputlist, TList* outputlist, const AliAODJet* jet, Double_t radius, Double_t& sumPt);
@@ -287,7 +287,7 @@ class AliAnalysisTaskFragmentationFunction : public AliAnalysisTaskSE {
   void     FillBckgHistos(Int_t type, TList* inputtracklist, TList* inputjetlist, AliAODJet* jet, 
                          AliFragFuncHistos* ffbckghistocuts,AliFragFuncQATrackHistos* qabckghistos,TH1F* fh1Mult = 0); 
  
-  Double_t GetMCStrangenessFactor(const Double_t& pt);
+  Double_t GetMCStrangenessFactor(const Double_t pt);
   Double_t GetMCStrangenessFactorCMS(AliAODMCParticle* daughter);
 
   void FillJetShape(AliAODJet* jet, TList* list,  TProfile* hProNtracksLeadingJet, TProfile** hProDelRPtSum, TProfile* hProDelR80pcPt=0, Double_t dPhiUE=0, Double_t normUE = 0, Bool_t scaleStrangeness = kFALSE);
@@ -332,6 +332,7 @@ class AliAnalysisTaskFragmentationFunction : public AliAnalysisTaskSE {
   UInt_t  fEvtSelectionMask;    // trigger class selection
   Int_t   fEventClass;          // centrality class selection
   Float_t fMaxVertexZ;          // maximum abs(z) position of primiary vertex [cm]
+  Bool_t  fRejectPileup;        // SPD pileup rejection
 
   // track cuts
   Float_t fTrackPtCut;    // track transverse momentum cut
index 32e5d56..40e5654 100644 (file)
@@ -66,6 +66,7 @@ ClassImp(AliAnalysisTaskJetProperties)
     ,fBranchJets("jets")
     ,fTrackType(kTrackAOD)
     ,fJetRejectType(0)
+    ,fRejectPileup(1)
     ,fUseAODInputJets(kTRUE)
     ,fFilterMask(0)
     ,fUsePhysicsSelection(kTRUE)
@@ -174,6 +175,7 @@ AliAnalysisTaskJetProperties::AliAnalysisTaskJetProperties(const char *name)
   ,fBranchJets("jets")
   ,fTrackType(kTrackAOD)
   ,fJetRejectType(0)
+  ,fRejectPileup(1)
   ,fUseAODInputJets(kTRUE)
   ,fFilterMask(0)
   ,fUsePhysicsSelection(kTRUE)
@@ -346,7 +348,7 @@ void AliAnalysisTaskJetProperties::UserCreateOutputObjects()
   fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected");
   fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected");
   fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected");
-  fh1EvtSelection->GetXaxis()->SetBinLabel(6,"vertex type: rejected");
+  fh1EvtSelection->GetXaxis()->SetBinLabel(6,"pile up: rejected");
   
   
   fh1VertexNContributors     = new TH1F("fh1VertexNContributors", "Vertex N contributors", 2500,-.5, 2499.5);
@@ -852,6 +854,14 @@ void AliAnalysisTaskJetProperties::UserExec(Option_t *)
     return;
   }
   
+  if(fRejectPileup && AliAnalysisHelperJetTasks::IsPileUp()){
+    if (fDebug > 2) Printf("%s:%d SPD pileup : event REJECTED...",(char*)__FILE__,__LINE__);
+    fh1EvtSelection->Fill(6.);
+    PostData(1, fCommonHistList);
+    return; 
+  }//pile up rejection
+  
+
   fh1VertexZ->Fill(primVtx->GetZ());
   
   if(TMath::Abs(primVtx->GetZ())>fMaxVertexZ){
index 40d97bd..0f5f00c 100644 (file)
@@ -44,6 +44,7 @@ class AliAnalysisTaskJetProperties : public AliAnalysisTaskSE {
   {fMaxVertexZ = VtxZ; fNContributors = nContributors;}
   virtual void   SetTrackCuts(Float_t trackPt = 0.15, Float_t trackEtaMin = -0.9, Float_t trackEtaMax = 0.9)
   {fTrackPtCut = trackPt; fTrackEtaMin = trackEtaMin; fTrackEtaMax = trackEtaMax;}
+  virtual void   SetPileupRejection(Bool_t IsPileupReject){fRejectPileup = IsPileupReject;}
   virtual void   SetJetCuts(Float_t jetPt = 5., Float_t jetEtaMin = -0.5, Float_t jetEtaMax = 0.5) 
   {fJetPtCut = jetPt; fJetEtaMin = jetEtaMin; fJetEtaMax = jetEtaMax;}
   virtual void   SetJetRejectType(Int_t i){fJetRejectType = i;}
@@ -53,7 +54,7 @@ class AliAnalysisTaskJetProperties : public AliAnalysisTaskSE {
   
   enum {kTrackUndef=0, kTrackAOD, kTrackKine,kTrackAODMC};//for track selection
   enum {kNoReject=0, kReject1Track};//for jet rejection
-  
+  enum {kRejectPileup=1};//for pileup rejection
  protected:
   Int_t           GetListOfJetTracks(TList* l, const AliAODJet* j);
   Int_t           GetListOfJets(TList* list);
@@ -75,6 +76,7 @@ class AliAnalysisTaskJetProperties : public AliAnalysisTaskSE {
   TString fBranchJets;          // branch name for reconstructed jets
   Int_t   fTrackType;           // type of generated tracks
   Int_t   fJetRejectType;       // type of jets rejected
+  Bool_t  fRejectPileup;        // pileup rejection
   Bool_t  fUseAODInputJets;     // take jets from in/output - only relevant if AOD event both in input AND output and we want to use output
   UInt_t  fFilterMask;         // filter bit for selected tracks
   Bool_t  fUsePhysicsSelection; // switch for event selection
index 3973295..885c299 100644 (file)
@@ -81,6 +81,10 @@ AliAnalysisTaskJetProperties *AddTaskJetProperties(Char_t* bJet="clustersAOD",
   else if(JetBranch.Contains("AOD"))   task->SetTrackType(AliAnalysisTaskJetProperties::kTrackAOD);
   else task->SetTrackType(AliAnalysisTaskJetProperties::kTrackUndef);//undefined track type
   
+  //setting pileup rejection
+  Bool_t IsPileupReject = kTRUE;//=kTRUE if you want to reject pileup, =kFALSE for no pileup rejection
+  task->SetPileupRejection(IsPileupReject);
+  
   //setting the jet rejection
   //two types implemented:kNoReject and kReject1Trk
   TString contName="";