hist switch, leading jet option

diff --git a/PWGGA/EMCALJetTasks/AliAnalysisTaskRho.cxx b/PWGGA/EMCALJetTasks/AliAnalysisTaskRho.cxx
index 44fd45dfc1bd470171e51b1338133ba08784fb72..54e564ed9cbb9eab81f3afffd31905f294ec4f76 100644 (file)
--- a/PWGGA/EMCALJetTasks/AliAnalysisTaskRho.cxx
+++ b/PWGGA/EMCALJetTasks/AliAnalysisTaskRho.cxx
@@ -4,12 +4,12 @@
 //
 // Authors: R.Reed, S.Aiola
 
-
 #include <TList.h>
 #include <TH1F.h>
 #include <TH2F.h>
 #include <TClonesArray.h>
 #include <TF1.h>
+#include <TMath.h>
 
 #include "AliLog.h"
 #include "AliAnalysisManager.h"
@@ -34,7 +34,9 @@ AliAnalysisTaskRho::AliAnalysisTaskRho() :
   fEtaMin(0),
   fEtaMax(0),
   fAreaCut(0),
+  fNExclLeadJets(0),
   fScaleFunction(0),
+  fCreateHisto(kFALSE),
   fOutputList(0),
   fHistCentrality(0),
   fHistJetPt(0),
@@ -51,8 +53,7 @@ AliAnalysisTaskRho::AliAnalysisTaskRho() :
   fHistJetPtvsNtrack(0),
   fHistJetAreavsNtrack(0),
   fHistNjetvsNtrack(0),
-  fRhoScaled(0),
-  fNewRhoFunction(0)
+  fRhoScaled(0)
 {
   // Constructor
 }
@@ -68,7 +69,9 @@ AliAnalysisTaskRho::AliAnalysisTaskRho(const char *name) :
   fEtaMin(0),
   fEtaMax(0),
   fAreaCut(0),
+  fNExclLeadJets(0),
   fScaleFunction(0),
+  fCreateHisto(kFALSE),
   fOutputList(0),
   fHistCentrality(0),
   fHistJetPt(0),
@@ -85,8 +88,7 @@ AliAnalysisTaskRho::AliAnalysisTaskRho(const char *name) :
   fHistJetPtvsNtrack(0),
   fHistJetAreavsNtrack(0),
   fHistNjetvsNtrack(0),
-  fRhoScaled(0),
-  fNewRhoFunction(0)
+  fRhoScaled(0)
 {
   // Constructor
 
@@ -100,12 +102,17 @@ void AliAnalysisTaskRho::UserCreateOutputObjects()
 
   fRhoScaledName = fRhoName;
   fRhoScaledName += "_Scaled";
-  fRhoScaled = new TParameter<Double_t>(fRhoScaledName, 0);
+  fRhoScaled = new TParameter<Double_t>(fRhoScaledName, 0);  
 
   OpenFile(1);
   fOutputList = new TList();
   fOutputList->SetOwner();
 
+  if (!fCreateHisto) {
+    PostData(1, fOutputList);
+    return;
+  }
+
   fHistCentrality             = new TH1F("Centrality",            "Centrality",            101, -1,  100);
   fHistRhovsCent              = new TH2F("RhovsCent",             "RhovsCent",             101, -1,  100,   500,  0,   500);
   fHistDeltaRhovsCent         = new TH2F("DeltaRhovsCent",        "DetlaRhovsCent",        101, -1,  100,   500, -250, 250);
@@ -123,8 +130,6 @@ void AliAnalysisTaskRho::UserCreateOutputObjects()
 
   fHistJetPt                  = new TH1F("JetPt",                  "Jet Pt",               100,   0,    250);
   fHistJetArea                = new TH1F("JetArea",                "Jet Area",             100, 0.0,    1.0);
-
-  fNewRhoFunction             = new TF1("rfunc","pol4", -1, 100);
   
   fOutputList->Add(fHistCentrality);
   fOutputList->Add(fHistRhovsCent);
@@ -142,8 +147,6 @@ void AliAnalysisTaskRho::UserCreateOutputObjects()
   fOutputList->Add(fHistJetPtvsNtrack);
   fOutputList->Add(fHistJetAreavsNtrack);
   fOutputList->Add(fHistNjetvsNtrack);
-
-  //fOutputList->Add(fNewRhoFunction);
   
   PostData(1, fOutputList);
   
@@ -158,6 +161,7 @@ Double_t AliAnalysisTaskRho::GetScaleFactor(Double_t cent)
   return scale;
 }
 
+
 //________________________________________________________________________
 void AliAnalysisTaskRho::UserExec(Option_t *) 
 {
@@ -165,11 +169,16 @@ void AliAnalysisTaskRho::UserExec(Option_t *)
   
   AliAnalysisTaskRhoBase::UserExec("");
 
-  // add rho to event if not yet there
+  fRho->SetVal(-1);
+
+   // add rho to event if not yet there
   if (!(InputEvent()->FindListObject(fRhoScaledName))) {
-    new(fRhoScaled) TParameter<Double_t>(fRhoScaledName, 0);
+    new(fRhoScaled) TParameter<Double_t>(fRhoScaledName, -1);
     InputEvent()->AddObject(fRhoScaled);
   }
+  else {
+    fRhoScaled->SetVal(-1);
+  }
 
   // optimization in case autobranch loading is off
   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
@@ -187,22 +196,62 @@ void AliAnalysisTaskRho::UserExec(Option_t *)
     
   jets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
   if (!jets) {
-    AliError(Form("Pointer to jets %s == 0", fTracksName.Data() ));
+    AliError(Form("Pointer to jets %s == 0", fJetsName.Data() ));
     return;
   }
-
-  fHistCentrality->Fill(fCent);
+  
+  if (fCreateHisto)
+    fHistCentrality->Fill(fCent);
 
   const Int_t Ntracks = tracks->GetEntries();
   const Int_t Njets = jets->GetEntries();
+
+  Int_t maxJetIds[] = {-1, -1};
+  Int_t maxJetPts[] = {0, 0};
+  if (fNExclLeadJets > 0) {
+    for (Int_t ij = 0; ij < Njets; ij++) {
+      
+      AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ij));
+      
+      if (!jet) {
+       AliError(Form("Could not receive jet %d", ij));
+       continue;
+      } 
+      
+      if (jet->Pt() > maxJetPts[0]) {
+       maxJetPts[1] = maxJetPts[0];
+       maxJetIds[1] = maxJetPts[0];
+       maxJetPts[0] = jet->Pt();
+       maxJetIds[0] = ij;
+      }
+      
+      if (jet->Pt() > maxJetPts[1]) {
+       maxJetPts[1] = jet->Pt();
+       maxJetIds[1] = ij;
+      }
+    }
+    if (fNExclLeadJets < 2) {
+      maxJetIds[1] = -1;
+      maxJetPts[1] = -1;
+    }
+  }
+
+  static Double_t rhovec[999];
   Int_t NjetAcc = 0;
 
-  vector<Double_t> rhovec;
+  // push all jets within selected acceptance into stack
   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
-    //push all jets within selected acceptance into stack
+    
+    // exlcuding lead jets
+    if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
+      continue;
+
     AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(iJets));
+
     if (!jet)
       continue; 
+
+    // applying some other cuts
     if (jet->Area() < fAreaCut)
       continue;
     if ((jet->Phi() < fPhiMin) || (jet->Phi() > fPhiMax))
@@ -211,100 +260,59 @@ void AliAnalysisTaskRho::UserExec(Option_t *)
       continue;
     if (jet->Area() == 0)
       continue;
+
+    rhovec[NjetAcc] = jet->Pt() / jet->Area();
+
     NjetAcc++;
-    rhovec.push_back(jet->Pt() / jet->Area());
-    fHistJetPt->Fill(jet->Pt());
-    fHistJetArea->Fill(jet->Area());
-    fHistJetPtvsCent->Fill(fCent, jet->Pt());
-    fHistJetPtvsNtrack->Fill(Ntracks, jet->Pt());
-    fHistJetAreavsCent->Fill(fCent, jet->Area());
-    fHistJetAreavsNtrack->Fill(Ntracks, jet->Area());
+
+    if (fCreateHisto) {
+      // filling histograms
+      fHistJetPt->Fill(jet->Pt());
+      fHistJetArea->Fill(jet->Area());
+      fHistJetPtvsCent->Fill(fCent, jet->Pt());
+      fHistJetPtvsNtrack->Fill(Ntracks, jet->Pt());
+      fHistJetAreavsCent->Fill(fCent, jet->Area());
+      fHistJetAreavsNtrack->Fill(Ntracks, jet->Area());
+    }
   }
   
-  fHistNjetvsCent->Fill(fCent,NjetAcc);
-  fHistNjetvsNtrack->Fill(Ntracks,NjetAcc);
+  if (fCreateHisto) {
+    fHistNjetvsCent->Fill(fCent, NjetAcc);
+    fHistNjetvsNtrack->Fill(Ntracks, NjetAcc);
+  }
 
   Double_t scale = GetScaleFactor(fCent);
   Double_t rhochem = GetRhoFactor(fCent);
 
   Double_t rho0 = -1;
   
-  if (rhovec.size() > 0){
+  if (NjetAcc > 0){
     //find median value
-    Sort(rhovec);
-    rho0 = GetMedian(rhovec,0);
-    fHistRhovsCent->Fill(fCent, rho0);
-    fHistDeltaRhovsCent->Fill(fCent, rho0 - rhochem);
-    fHistDeltaRhoScalevsCent->Fill(fCent, rho0 * scale - rhochem);
-    fHistRhovsNtrack->Fill(Ntracks, rho0);
-    fHistDeltaRhovsNtrack->Fill(Ntracks, rho0 - rhochem);
-    fHistDeltaRhoScalevsNtrack->Fill(Ntracks, rho0 * scale - rhochem);
-  }
-
-  fRho->SetVal(rho0);
+    rho0 = TMath::Median(NjetAcc, rhovec);
 
-  Double_t rho_scaled = rho0 * scale;
+    Double_t rhoScaled = rho0 * scale;
 
-  fRhoScaled->SetVal(rho_scaled);
-  
-  PostData(1, fOutputList);
-}      
+    fRho->SetVal(rho0);
+    fRhoScaled->SetVal(rhoScaled);
 
-//________________________________________________________________________
-void AliAnalysisTaskRho::Sort(vector<Double_t>& v)
-{
-  vector<Double_t> temp;
-  temp.push_back(v[0]);
-  for (UInt_t i = 1; i < v.size(); i++) {
-    Bool_t insert = kFALSE;
-    for (vector<Double_t>::iterator j = temp.begin(); j < temp.end(); j++){
-      if (v[i]> * j) {
-       temp.insert(j, v[i]);
-       insert = kTRUE;
-       j = temp.end();
-      }
+    if (fCreateHisto) {
+      // filling other histograms
+      fHistRhovsCent->Fill(fCent, rho0);
+      fHistDeltaRhovsCent->Fill(fCent, rho0 - rhochem);
+      fHistDeltaRhoScalevsCent->Fill(fCent, rhoScaled - rhochem);
+      fHistRhovsNtrack->Fill(Ntracks, rho0);
+      fHistDeltaRhovsNtrack->Fill(Ntracks, rho0 - rhochem);
+      fHistDeltaRhoScalevsNtrack->Fill(Ntracks, rhoScaled - rhochem);
     }
-    if (!insert)
-      temp.insert(temp.end(), v[i]);
   }
-  v = temp;
-  return;
-}
 
-//________________________________________________________________________
-Double_t AliAnalysisTaskRho::GetMedian(vector<Double_t> v, Int_t c)
-{
-  if (v.size() == 0)
-    return -1;
-  if ((v.size() < 2) && (c == 1))
-    return -1;
-  if ((v.size() < 3) && ( c== 2))
-    return -1;
-  
-  if (c == 1)
-    v.erase(v.begin());
-  
-  if (c == 2) {
-    v.erase(v.begin());
-    v.erase(v.begin());
-  }
+  if (fCreateHisto)
+    PostData(1, fOutputList);
+}      
 
-  Float_t middle = (Float_t)v.size() / 2.0;
-  if (middle != ceil(middle)) {
-    //odd number
-    return v[floor(middle)];
-  }
-  else{
-    //even
-    return (v[middle] + v[middle-1]) / 2.0;
-  } 
-}
 
 //________________________________________________________________________
 void AliAnalysisTaskRho::Terminate(Option_t *) 
 {
-  /*
-  fHistRhovsCent->Fit(fNewRhoFunction, "NO");
-  PostData(1, fOutputList);
-  */
+
 }

diff --git a/PWGGA/EMCALJetTasks/AliAnalysisTaskRho.h b/PWGGA/EMCALJetTasks/AliAnalysisTaskRho.h
index d483c327be0a7c1e66225f39542509b1d14e69e0..db14abece8b3b91b6c5376f67fa2e82347904a81 100644 (file)
--- a/PWGGA/EMCALJetTasks/AliAnalysisTaskRho.h
+++ b/PWGGA/EMCALJetTasks/AliAnalysisTaskRho.h
@@ -31,11 +31,11 @@ class AliAnalysisTaskRho : public AliAnalysisTaskRhoBase {
   void                   SetJetsName(const char *n)                            { fJetsName      = n    ; }
   void                   SetScaleFunction(TF1* sf)                             { fScaleFunction = sf   ; }
   void                   SetTracksName(const char *n)                          { fTracksName    = n    ; }
+  void                   SetExcludeLeadJets(UInt_t n)                          { fNExclLeadJets = n    ; }
+  void                   SetCreateHistograms(Bool_t h = kTRUE)                 { fCreateHisto   = h    ; }
   
  protected:
-  virtual void           Sort(vector<Double_t>& v)            ;
-  virtual Double_t       GetMedian(vector<Double_t> v, int c) ;
-  virtual Double_t       GetScaleFactor(Double_t cent)        ;
+  virtual Double_t       GetScaleFactor(Double_t cent);
 
   TString                fTracksName;                    // name of track collection
   TString                fJetsName;                      // name of jet collection
@@ -46,7 +46,9 @@ class AliAnalysisTaskRho : public AliAnalysisTaskRhoBase {
   Double_t               fEtaMin;                        // minimum eta
   Double_t               fEtaMax;                        // maximum eta
   Double_t               fAreaCut;                       // cut on jet area
+  UInt_t                 fNExclLeadJets;                 // number of leading jets to be excluded from the median calculation
   TF1                   *fScaleFunction;                 // pre-computed scale factor as a function of centrality
+  Bool_t                 fCreateHisto;                   // whether or not create histograms
   TList                 *fOutputList;                    //!output list
   TH1F                  *fHistCentrality;                //!centrality distribution
   TH1F                  *fHistJetPt;                     //!jet pt distribution
@@ -65,7 +67,6 @@ class AliAnalysisTaskRho : public AliAnalysisTaskRhoBase {
   TH2F                  *fHistJetAreavsNtrack;           //!jet area vs. no. of tracks
   TH2F                  *fHistNjetvsNtrack;              //!no. of jets vs. no. of tracks
   TParameter<Double_t>  *fRhoScaled;                     //!per event scaled rho
-  TF1                   *fNewRhoFunction;                //!new rho as a function of centrality
 
   AliAnalysisTaskRho(const AliAnalysisTaskRho&);             // not implemented
   AliAnalysisTaskRho& operator=(const AliAnalysisTaskRho&);  // not implemented

diff --git a/PWGGA/EMCALJetTasks/macros/AddTaskRho.C b/PWGGA/EMCALJetTasks/macros/AddTaskRho.C
index d1d4036af24cec0eda8bb4a07e55054a3b97ad1d..0eb43493e0e0666aa88ad21e7ef03ad28038952e 100644 (file)
--- a/PWGGA/EMCALJetTasks/macros/AddTaskRho.C
+++ b/PWGGA/EMCALJetTasks/macros/AddTaskRho.C
@@ -9,7 +9,8 @@ AliAnalysisTaskRho* AddTaskRho(
    const Double_t maxPhi      = 2 * TMath::Pi(),
    const Double_t minEta      = -0.3,
    const Double_t maxEta      = 0.3,
-   const Double_t minArea     = 0.0
+   const Double_t minArea     = 0.0,
+   const UInt_t   exclJets    = 0
 )
 {  
   // Get the pointer to the existing analysis manager via the static access method.
@@ -41,6 +42,7 @@ AliAnalysisTaskRho* AddTaskRho(
   rhotask->SetJetPhi(minPhi,maxPhi);
   rhotask->SetJetEta(minEta,maxEta);
   rhotask->SetAreaCut(minArea);
+  rhotask->SetExcludeLeadJets(exclJets);
 
   //-------------------------------------------------------
   // Final settings, pass to manager and set the containers