rho for sparse systems from megan

diff --git a/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg b/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg
index 4944378669ee313f61f70fa2baeb888171c6f1e7..c6c13393089f04c0779635f55bed48fe40ff92e7 100644 (file)
--- a/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg
+++ b/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg
@@ -33,6 +33,7 @@ set ( SRCS
  EMCALJetTasks/AliAnalysisTaskRhoAverage.cxx
  EMCALJetTasks/AliAnalysisTaskRhoBase.cxx
  EMCALJetTasks/AliAnalysisTaskRhoFlow.cxx
+ EMCALJetTasks/AliAnalysisTaskRhoSparse.cxx
  EMCALJetTasks/AliAnalysisTaskScale.cxx
  EMCALJetTasks/AliEmcalJet.cxx
  EMCALJetTasks/AliEmcalJetTask.cxx

diff --git a/PWGJE/EMCALJetTasks/AliAnalysisTaskRhoSparse.cxx b/PWGJE/EMCALJetTasks/AliAnalysisTaskRhoSparse.cxx
new file mode 100644 (file)
index 0000000..41ef129
--- /dev/null
+++ b/PWGJE/EMCALJetTasks/AliAnalysisTaskRhoSparse.cxx
@@ -0,0 +1,157 @@
+// $Id: AliAnalysisTaskRho.cxx 58408 2012-09-03 07:00:58Z loizides $
+//
+// Calculation of rho from a collection of jets.
+// If scale function is given the scaled rho will be exported
+// with the name as "fRhoName".apppend("_Scaled").
+//
+// Authors: R.Reed, S.Aiola, M.Connors
+
+#include "AliAnalysisTaskRhoSparse.h"
+
+#include <TClonesArray.h>
+#include <TMath.h>
+
+#include "AliAnalysisManager.h"
+#include "AliEmcalJet.h"
+#include "AliLog.h"
+#include "AliRhoParameter.h"
+
+ClassImp(AliAnalysisTaskRhoSparse)
+
+//________________________________________________________________________
+AliAnalysisTaskRhoSparse::AliAnalysisTaskRhoSparse() : 
+  AliAnalysisTaskRhoBase("AliAnalysisTaskRhoSparse"),
+  fHistOccCorrvsCent(0),
+  fNExclLeadJets(0),
+  fRhoCMS(0)
+{
+  // Constructor.
+}
+
+//________________________________________________________________________
+AliAnalysisTaskRhoSparse::AliAnalysisTaskRhoSparse(const char *name, Bool_t histo) :
+  AliAnalysisTaskRhoBase(name, histo),
+  fHistOccCorrvsCent(0),
+  fNExclLeadJets(0),
+  fRhoCMS(0)
+{
+  // Constructor.
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskRhoSparse::UserCreateOutputObjects()
+{
+  if (!fCreateHisto)
+    return;
+
+  AliAnalysisTaskRhoBase::UserCreateOutputObjects();
+  
+  fHistOccCorrvsCent = new TH2F("OccCorrvsCent", "OccCorrvsCent", 101, -1, 100, 2000, 0 , 2);
+  fOutput->Add(fHistOccCorrvsCent);
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskRhoSparse::Run() 
+{
+  // Run the analysis.
+
+  fRho->SetVal(0);
+  if (fRhoScaled)
+    fRhoScaled->SetVal(0);
+
+  if (!fJets)
+    return kFALSE;
+
+  const Int_t Njets   = fJets->GetEntries();
+
+  Int_t maxJetIds[]   = {-1, -1};
+  Float_t maxJetPts[] = { 0,  0};
+
+  if (fNExclLeadJets > 0) {
+    for (Int_t ij = 0; ij < Njets; ++ij) {
+      AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
+      if (!jet) {
+       AliError(Form("%s: Could not receive jet %d", GetName(), ij));
+       continue;
+      } 
+
+      if (!AcceptJet(jet))
+        continue;
+
+      if (jet->Pt() > maxJetPts[0]) {
+       maxJetPts[1] = maxJetPts[0];
+       maxJetIds[1] = maxJetIds[0];
+       maxJetPts[0] = jet->Pt();
+       maxJetIds[0] = ij;
+      } else if (jet->Pt() > maxJetPts[1]) {
+       maxJetPts[1] = jet->Pt();
+       maxJetIds[1] = ij;
+      }
+    }
+    if (fNExclLeadJets < 2) {
+      maxJetIds[1] = -1;
+      maxJetPts[1] = 0;
+    }
+  }
+
+  static Double_t rhovec[999];
+  Int_t NjetAcc = 0;
+  Double_t TotaljetArea=0;
+  Double_t TotaljetAreaPhys=0;
+
+  // push all jets within selected acceptance into stack
+  for (Int_t iJets = 0; iJets < Njets; ++iJets) {
+
+    // exlcuding lead jets
+    if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
+      continue;
+
+    AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
+    if (!jet) {
+      AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
+      continue;
+    } 
+
+    //cout << "jetpt: " << jet->Pt() << " jetArea: " << jet->Area() <<endl;
+
+    if (!AcceptJet(jet))
+      continue;
+
+
+    if(jet->Pt()>0.01) rhovec[NjetAcc] = jet->Pt() / jet->Area();
+    //cout << "ACCEPTED: jetpt: " << jet->Pt() << " jetArea: " << jet->Area() << " jetRho: " << rhovec[NjetAcc] <<endl;
+    if(jet->Pt()>0.01) TotaljetAreaPhys+=jet->Area();
+    TotaljetArea+=jet->Area();
+    ++NjetAcc;
+  }
+
+  const Double_t TpcMaxPhi = TMath::Pi()*2.;
+
+  const Double_t TpcArea = TpcMaxPhi * 2.*(0.7);
+  Double_t OccCorr=0.0;
+  //cout << "Area Physical: " << TotaljetAreaPhys << " total: " << TotaljetArea <<endl;
+  //if(TotaljetArea>0) OccCorr=TotaljetAreaPhys/TotaljetArea;
+  if(TpcArea>0) OccCorr=TotaljetAreaPhys/TpcArea;
+ 
+  fHistOccCorrvsCent->Fill(fCent, OccCorr);
+
+
+  if (NjetAcc > 0) {
+    //find median value
+    Double_t rho = TMath::Median(NjetAcc, rhovec);
+
+    if(fRhoCMS){
+      rho = rho * OccCorr;
+    }
+
+    fRho->SetVal(rho);
+
+
+    if (fRhoScaled) {
+      Double_t rhoScaled = rho * GetScaleFactor(fCent);
+      fRhoScaled->SetVal(rhoScaled);
+    }
+  }
+
+  return kTRUE;
+} 

diff --git a/PWGJE/EMCALJetTasks/AliAnalysisTaskRhoSparse.h b/PWGJE/EMCALJetTasks/AliAnalysisTaskRhoSparse.h
new file mode 100644 (file)
index 0000000..3cb4fda
--- /dev/null
+++ b/PWGJE/EMCALJetTasks/AliAnalysisTaskRhoSparse.h
@@ -0,0 +1,34 @@
+#ifndef ALIANALYSISTASKRHOSPARSE_H
+#define ALIANALYSISTASKRHOSPARSE_H
+
+// $Id: AliAnalysisTaskRho.h 58408 2012-09-03 07:00:58Z loizides $
+
+#include "AliAnalysisTaskRhoBase.h"
+
+class AliAnalysisTaskRhoSparse : public AliAnalysisTaskRhoBase {
+
+ public:
+  AliAnalysisTaskRhoSparse();
+  AliAnalysisTaskRhoSparse(const char *name, Bool_t histo=kFALSE);
+  virtual ~AliAnalysisTaskRhoSparse() {}
+
+  void             UserCreateOutputObjects();
+  void             SetExcludeLeadJets(UInt_t n)    { fNExclLeadJets = n    ; }
+  void             SetRhoCMS(Bool_t cms)           { fRhoCMS = cms ; }
+
+
+ protected:
+  Bool_t           Run();
+  TH2F             *fHistOccCorrvsCent;             //!occupancy correction vs. centrality
+
+  UInt_t           fNExclLeadJets;                 // number of leading jets to be excluded from the median calculation
+
+  Bool_t           fRhoCMS;
+
+
+  AliAnalysisTaskRhoSparse(const AliAnalysisTaskRhoSparse&);             // not implemented
+  AliAnalysisTaskRhoSparse& operator=(const AliAnalysisTaskRhoSparse&);  // not implemented
+  
+  ClassDef(AliAnalysisTaskRhoSparse, 1); // Rho task
+};
+#endif

diff --git a/PWGJE/EMCALJetTasks/macros/AddTaskRhoSparse.C b/PWGJE/EMCALJetTasks/macros/AddTaskRhoSparse.C
new file mode 100644 (file)
index 0000000..7bf3daf
--- /dev/null
+++ b/PWGJE/EMCALJetTasks/macros/AddTaskRhoSparse.C
@@ -0,0 +1,81 @@
+// $Id: AddTaskRhoSparse.C 58584 2012-09-13 10:37:42Z loizides $
+
+AliAnalysisTaskRhoSparse* AddTaskRhoSparse(
+   const char    *nJets       = "Jets",
+   const char    *nTracks     = "PicoTracks",
+   const char    *nClusters   = "CaloClusters",  
+   const char    *nRho        = "Rho",
+   Double_t       jetradius   = 0.2,
+   UInt_t         type        = AliAnalysisTaskEmcal::kTPC,
+   Double_t       jetareacut  = 0.01,
+   Double_t       jetptcut    = 0.0,
+   Double_t       emcareacut  = 0,
+   TF1           *sfunc       = 0,
+   const UInt_t   exclJets    = 2,
+   const Bool_t   histo       = kFALSE,
+   const char    *taskname    = "Rho",
+   const Bool_t   fRhoCMS      = kTRUE
+)
+{  
+  // Get the pointer to the existing analysis manager via the static access method.
+  //==============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr)
+  {
+    ::Error("AddTaskRhoSparse", "No analysis manager to connect to.");
+    return NULL;
+  }  
+  
+  // Check the analysis type using the event handlers connected to the analysis manager.
+  //==============================================================================
+  if (!mgr->GetInputEventHandler())
+  {
+    ::Error("AddTaskRhoSparse", "This task requires an input event handler");
+    return NULL;
+  }
+  
+  //-------------------------------------------------------
+  // Init the task and do settings
+  //-------------------------------------------------------
+
+  TString name(Form("%s_%s_", taskname, nJets));
+  if (type == AliAnalysisTaskEmcal::kTPC) 
+    name += "TPC";
+  else if (type == AliAnalysisTaskEmcal::kEMCAL) 
+    name += "EMCAL";
+  else if (type == AliAnalysisTaskEmcal::kUser) 
+    name += "USER";
+  AliAnalysisTaskRhoSparse *rhotask = new AliAnalysisTaskRhoSparse(name, histo);
+  rhotask->SetHistoBins(1000,-0.1,9.9);
+  rhotask->SetAnaType(type);
+  rhotask->SetScaleFunction(sfunc);
+  rhotask->SetJetsName(nJets);
+  rhotask->SetTracksName(nTracks);
+  rhotask->SetClusName(nClusters);
+  rhotask->SetRhoName(nRho);
+  rhotask->SetJetAreaCut(jetareacut);
+  rhotask->SetAreaEmcCut(emcareacut);
+  rhotask->SetJetPtCut(jetptcut);
+  rhotask->SetJetRadius(jetradius);
+  rhotask->SetExcludeLeadJets(exclJets);
+  rhotask->SetRhoCMS(fRhoCMS);
+
+  //-------------------------------------------------------
+  // Final settings, pass to manager and set the containers
+  //-------------------------------------------------------
+
+  mgr->AddTask(rhotask);
+
+  // Create containers for input/output
+  mgr->ConnectInput(rhotask, 0, mgr->GetCommonInputContainer());
+  if (histo) {
+    TString contname(name);
+    contname += "_histos";
+    AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(), 
+                                                             TList::Class(),AliAnalysisManager::kOutputContainer,
+                                                             Form("%s", AliAnalysisManager::GetCommonFileName()));
+    mgr->ConnectOutput(rhotask, 1, coutput1);
+  }
+
+  return rhotask;
+}

diff --git a/PWGJE/PWGJEEMCALJetTasksLinkDef.h b/PWGJE/PWGJEEMCALJetTasksLinkDef.h
index 16bba9044828f5e43934fb45dc6e6464a4f1d446..456afc602825f2f2a65c075b5613d34a0e32219a 100644 (file)
--- a/PWGJE/PWGJEEMCALJetTasksLinkDef.h
+++ b/PWGJE/PWGJEEMCALJetTasksLinkDef.h
@@ -11,6 +11,7 @@
 #pragma link C++ class AliAnalysisTaskRho+;
 #pragma link C++ class AliAnalysisTaskRhoFlow+;
 #pragma link C++ class AliAnalysisTaskRhoAverage+;
+#pragma link C++ class AliAnalysisTaskRhoSparse+;
 #pragma link C++ class AliAnalysisTaskDeltaPt+;
 #pragma link C++ class AliAnalysisTaskScale+;
 #pragma link C++ class AliEmcalJet+;