updates from pdsf
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 17 May 2012 02:51:56 +0000 (02:51 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 17 May 2012 02:51:56 +0000 (02:51 +0000)
PWGGA/EMCALJetTasks/AliAnalysisTaskRho.cxx
PWGGA/EMCALJetTasks/AliAnalysisTaskRho.h

index b455f47..5249add 100644 (file)
@@ -1,9 +1,15 @@
-#include <TChain.h>
+// $Id: $
+//
+// Calculation of rho 
+//
+// Authors: R.Reed, S.Aiola
+
+
 #include <TList.h>
 #include <TH1F.h>
 #include <TH2F.h>
 #include <TClonesArray.h>
-#include <TVector3.h>
+#include <TF1.h>
 
 #include "AliLog.h"
 #include "AliAnalysisManager.h"
@@ -18,90 +24,83 @@ ClassImp(AliAnalysisTaskRho)
 
 //________________________________________________________________________
 AliAnalysisTaskRho::AliAnalysisTaskRho() : 
-  AliAnalysisTaskSE(),
+  AliAnalysisTaskRhoBase(),
   fTracksName("tracks"),
-  fJetsName("jets"),
-  fRhosName("fArrRhos"),
-  fClustersName("clusters"),
+  fJetsName("KtJets"),
+  fClustersName("caloClusters"),
+  fRhoScaledName(""),
+  fPhiMin(0),
+  fPhiMax(0),
+  fEtaMin(0),
+  fEtaMax(0),
+  fAreaCut(0),
+  fScaleFunction(0),
   fOutputList(0),
   fHistCentrality(0),
   fHistJetPt(0),
-  fHistJetArea(0), 
+  fHistJetArea(0),
   fHistRhovsCent(0),
-  fHistDeltaRhovsCent(0),  
-  fHistDeltaRhoScalevsCent(0),  
+  fHistDeltaRhovsCent(0),
+  fHistDeltaRhoScalevsCent(0),
   fHistJetPtvsCent(0),
-  fHistJetAreavsCent(0),  
+  fHistJetAreavsCent(0),
   fHistNjetvsCent(0),
   fHistRhovsNtrack(0),
-  fHistDeltaRhovsNtrack(0),   
-  fHistDeltaRhoScalevsNtrack(0),   
+  fHistDeltaRhovsNtrack(0),
+  fHistDeltaRhoScalevsNtrack(0),
   fHistJetPtvsNtrack(0),
   fHistJetAreavsNtrack(0),
   fHistNjetvsNtrack(0),
-  fArrRhos(0),
-  fPhiMin(0),
-  fPhiMax(0),
-  fEtaMin(0),
-  fEtaMax(0),
-  fAreaCut(0),
-  fCswitch(0)
+  fRhoScaled(0),
+  fNewRhoFunction(0)
 {
   // Constructor
-
-  DefineInput(0, TChain::Class());
-  DefineOutput(1, TList::Class());
 }
 //________________________________________________________________________
 AliAnalysisTaskRho::AliAnalysisTaskRho(const char *name) :
-  AliAnalysisTaskSE(name),
+  AliAnalysisTaskRhoBase(name),
   fTracksName("tracks"),
-  fJetsName("jets"),
-  fRhosName("fArrRhos"),
-  fClustersName("clusters"),
+  fJetsName("KtJets"),
+  fClustersName("caloClusters"),
+  fRhoScaledName(""),
+  fPhiMin(0),
+  fPhiMax(0),
+  fEtaMin(0),
+  fEtaMax(0),
+  fAreaCut(0),
+  fScaleFunction(0),
   fOutputList(0),
   fHistCentrality(0),
   fHistJetPt(0),
-  fHistJetArea(0), 
+  fHistJetArea(0),
   fHistRhovsCent(0),
-  fHistDeltaRhovsCent(0),  
-  fHistDeltaRhoScalevsCent(0),  
+  fHistDeltaRhovsCent(0),
+  fHistDeltaRhoScalevsCent(0),
   fHistJetPtvsCent(0),
-  fHistJetAreavsCent(0),  
+  fHistJetAreavsCent(0),
   fHistNjetvsCent(0),
   fHistRhovsNtrack(0),
-  fHistDeltaRhovsNtrack(0),   
-  fHistDeltaRhoScalevsNtrack(0),   
+  fHistDeltaRhovsNtrack(0),
+  fHistDeltaRhoScalevsNtrack(0),
   fHistJetPtvsNtrack(0),
   fHistJetAreavsNtrack(0),
   fHistNjetvsNtrack(0),
-  fArrRhos(0),
-  fPhiMin(0),
-  fPhiMax(0),
-  fEtaMin(0),
-  fEtaMax(0),
-  fAreaCut(0),
-  fCswitch(0)
-
+  fRhoScaled(0),
+  fNewRhoFunction(0)
 {
   // Constructor
 
-  DefineInput(0, TChain::Class());
   DefineOutput(1, TList::Class());
 }
 
 //________________________________________________________________________
 void AliAnalysisTaskRho::UserCreateOutputObjects()
 {
+  AliAnalysisTaskRhoBase::UserCreateOutputObjects();
 
-  AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
-  if (!handler) {
-    AliError("Input handler not available!");
-    return;
-  }
-
-  fArrRhos = new TClonesArray("TVector3");
-  fArrRhos->SetName(fRhosName);
+  fRhoScaledName = fRhoName;
+  fRhoScaledName += "_Scaled";
+  fRhoScaled = new TParameter<Double_t>(fRhoScaledName, 0);
 
   OpenFile(1);
   fOutputList = new TList();
@@ -124,6 +123,8 @@ 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);
@@ -141,6 +142,8 @@ void AliAnalysisTaskRho::UserCreateOutputObjects()
   fOutputList->Add(fHistJetPtvsNtrack);
   fOutputList->Add(fHistJetAreavsNtrack);
   fOutputList->Add(fHistNjetvsNtrack);
+
+  fOutputList->Add(fNewRhoFunction);
   
   PostData(1, fOutputList);
   
@@ -149,56 +152,40 @@ void AliAnalysisTaskRho::UserCreateOutputObjects()
 //________________________________________________________________________
 Double_t AliAnalysisTaskRho::GetScaleFactor(Double_t cent)
 {
-  Double_t scale = 0.000066*cent*cent-0.0015*cent+1.5;
+  Double_t scale = 1;
+  if (fScaleFunction)
+    scale = fScaleFunction->Eval(cent);
   return scale;
 }
 
 //________________________________________________________________________
-Double_t AliAnalysisTaskRho::GetRhoFactor(Double_t cent)
-{
-  Double_t RhoChEmArray[100] = {201.32,  191.28,  183.72,  176.70,  170.10,  163.63,  157.76,  152.13,  146.62,  141.28,  136.41,  131.69,  126.80,  121.95,  117.85,  113.44,  109.19,  105.35,  101.24,  97.23,  93.49,  90.15,  86.77,  83.26,  79.92,  76.77,  73.75,  70.92,  67.93,  65.09,  62.39,  59.57,  57.26,  54.65,  52.24,  50.03,  47.70,  45.34,  43.29,  41.32,  39.35,  37.35,  35.46,  33.64,  32.15,  30.43,  28.81,  27.28,  25.85,  24.36,  23.09,  21.83,  20.60,  19.42,  18.45,  17.44,  16.62,  15.40,  14.06,  13.31,  12.72,  11.91,  11.01,  10.54,  9.89,  9.23,  8.73,  8.15,  7.73,  7.17,  7.00,  6.55,  6.07,  5.78,  5.42,  5.29,  5.13,  4.69,  4.74,  4.36,  4.40,  4.21,  4.23,  3.77,  3.83,  4.23,  4.17,  4.53,  4.06,  3.63,  3.90,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00,  0.00};
-  Double_t rho = RhoChEmArray[(Int_t)floor(cent)];
-
-  return rho;
-}
-
-//________________________________________________________________________
 void AliAnalysisTaskRho::UserExec(Option_t *) 
 {
   // Main loop, called for each event.
+  
+  AliAnalysisTaskRhoBase::UserExec("");
 
-  // add rhos to event if not yet there
-  if (!(InputEvent()->FindListObject(fRhosName)))
-    InputEvent()->AddObject(fArrRhos);
+  // add rho to event if not yet there
+  if (!(InputEvent()->FindListObject(fRhoScaledName))) {
+    new(fRhoScaled) TParameter<Double_t>(fRhoScaledName, 0);
+    InputEvent()->AddObject(fRhoScaled);
+  }
 
   // optimization in case autobranch loading is off
   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
   if (fTracksName == "Tracks")
     am->LoadBranch("Tracks");
 
-  // get centrality 
-  Double_t fCent = -1; 
-  AliCentrality *centrality = InputEvent()->GetCentrality() ;
-  TList *l = InputEvent()->GetList();
-  if (centrality)
-    fCent = centrality->GetCentralityPercentile("V0M");
-  else
-    fCent=99; // probably pp data
-  if (fCent<0) {
-    AliError(Form("Centrality negative: %f", fCent));
-    return;
-  }
-
   TClonesArray *jets = 0;
   TClonesArray *tracks = 0;
 
-  tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
+  tracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
   if (!tracks) {
   AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
    return;
   }
     
-  jets = dynamic_cast<TClonesArray*>(l->FindObject(fJetsName));
+  jets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
   if (!jets) {
     AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
     return;
@@ -206,14 +193,11 @@ void AliAnalysisTaskRho::UserExec(Option_t *)
 
   fHistCentrality->Fill(fCent);
 
-  Double_t scale = GetScaleFactor(fCent);
-  Double_t rhochem = GetRhoFactor(fCent);
-
   const Int_t Ntracks = tracks->GetEntries();
   const Int_t Njets = jets->GetEntries();
   Int_t NjetAcc = 0;
 
-   vector<Double_t> rhovec;
+  vector<Double_t> rhovec;
   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
     //push all jets within selected acceptance into stack
     AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(iJets));
@@ -228,51 +212,45 @@ void AliAnalysisTaskRho::UserExec(Option_t *)
     if (jet->Area() == 0)
       continue;
     NjetAcc++;
-    rhovec.push_back(jet->Pt()/jet->Area());
+    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());
+    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);
+
+  Double_t scale = GetScaleFactor(fCent);
+  Double_t rhochem = GetRhoFactor(fCent);
+
   Double_t rho0 = -1;
   
-  if (rhovec.size()>0){
+  if (rhovec.size() > 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);
+    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);
   }
-  
-  //Scale rho for method 2, task does not know whether jet collection
-  //is charged only
 
-  Double_t rho2=scale*rho0;
-    
-  Int_t irx = 0;
-  TVector3 *ArrRhos = new ((*fArrRhos)[irx]) TVector3;
-  ArrRhos->SetXYZ(rho0, rho2, rhochem);
-  irx++;
+  fRho->SetVal(rho0);
+
+  Double_t rho_scaled = rho0 * scale;
+
+  fRhoScaled->SetVal(rho_scaled);
   
   PostData(1, fOutputList);
 }      
 
 //________________________________________________________________________
-void AliAnalysisTaskRho::Terminate(Option_t *) 
-{
-
-}
-
-//________________________________________________________________________
 void AliAnalysisTaskRho::Sort(vector<Double_t>& v)
 {
   vector<Double_t> temp;
@@ -321,3 +299,10 @@ Double_t AliAnalysisTaskRho::GetMedian(vector<Double_t> v, Int_t c)
     return (v[middle] + v[middle-1]) / 2.0;
   } 
 }
+
+//________________________________________________________________________
+void AliAnalysisTaskRho::Terminate(Option_t *) 
+{
+  fHistRhovsCent->Fit(fNewRhoFunction, "NO");
+  PostData(1, fOutputList);
+}
index 44f7d61..d483c32 100644 (file)
@@ -1,75 +1,75 @@
 #ifndef ALIANALYSISTASKRHO_cxx
 #define ALIANALYSISTASKRHO_cxx
 
+// $Id: $
+
 class TList;
 class TH1F;
 class TH2F;
 class TClonesArray;
 class TString;
+class TF1;
+
+#include <TParameter.h>
 
-#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisTaskRhoBase.h"
 
-class AliAnalysisTaskRho : public AliAnalysisTaskSE {
+class AliAnalysisTaskRho : public AliAnalysisTaskRhoBase {
 
  public:
   AliAnalysisTaskRho();
   AliAnalysisTaskRho(const char *name);
   virtual ~AliAnalysisTaskRho() {}
   
-  virtual void          UserCreateOutputObjects();
-  virtual void          UserExec(Option_t*);
-  virtual void          Terminate(Option_t*);
+  virtual void           UserCreateOutputObjects();
+  virtual void           UserExec(Option_t*);
+  virtual void           Terminate(Option_t*);
 
-  void                  SetTracksName(const char *n)                          { fTracksName  = n    ; }
-  void                  SetJetsName(const char *n)                            { fJetsName    = n    ; }
-  void                  SetRhosName(const char *rn)                           { fRhosName    = rn   ; }
-  void                  SetJetPhi(Double_t pmin = 1.8, Double_t pmax = 2.74)  { fPhiMin      = pmin ; fPhiMax = pmax  ; }
-  void                  SetJetEta(Double_t emin = -0.3, Double_t emax = 0.3)  { fEtaMin      = emin ; fEtaMax = emax  ; }
-  void                  SetAreaCut(Double_t a = 0.0)                          { fAreaCut     = a    ; }
+  void                   SetAreaCut(Double_t a = 0.0)                          { fAreaCut       = a    ; }
+  void                   SetJetEta(Double_t emin, Double_t emax)               { fEtaMin        = emin ; fEtaMax = emax  ; }
+  void                   SetJetPhi(Double_t pmin, Double_t pmax)               { fPhiMin        = pmin ; fPhiMax = pmax  ; }
+  void                   SetJetsName(const char *n)                            { fJetsName      = n    ; }
+  void                   SetScaleFunction(TF1* sf)                             { fScaleFunction = sf   ; }
+  void                   SetTracksName(const char *n)                          { fTracksName    = n    ; }
   
  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      GetRhoFactor(Double_t cent)          ;
+  virtual void           Sort(vector<Double_t>& v)            ;
+  virtual Double_t       GetMedian(vector<Double_t> v, int c) ;
+  virtual Double_t       GetScaleFactor(Double_t cent)        ;
 
- private: 
   TString                fTracksName;                    // name of track collection
   TString                fJetsName;                      // name of jet collection
-  TString                fRhosName;                      // name of Rho array output
   TString                fClustersName;                  // name of clusters collection
-
-  TList                 *fOutputList;                    //! Output list
-  TH1F                  *fHistCentrality;
-  TH1F                  *fHistJetPt;
-  TH1F                  *fHistJetArea; 
-  TH2F                  *fHistRhovsCent;
-  TH2F                  *fHistDeltaRhovsCent;  
-  TH2F                  *fHistDeltaRhoScalevsCent;  
-  TH2F                  *fHistJetPtvsCent;  
-  TH2F                  *fHistJetAreavsCent;  
-  TH2F                  *fHistNjetvsCent;
-
-  TH2F                 *fHistRhovsNtrack;
-  TH2F                 *fHistDeltaRhovsNtrack;   
-  TH2F                 *fHistDeltaRhoScalevsNtrack;   
-  TH2F                 *fHistJetPtvsNtrack;
-  TH2F                 *fHistJetAreavsNtrack;
-  TH2F                 *fHistNjetvsNtrack;
-
-  TClonesArray         *fArrRhos;
-
-  Double_t              fPhiMin;
-  Double_t              fPhiMax;
-  Double_t              fEtaMin;
-  Double_t              fEtaMax;
-  Double_t              fAreaCut;
-  Int_t                 fCswitch;
+  TString                fRhoScaledName;                 // name of scaled rho object
+  Double_t               fPhiMin;                        // minimum phi
+  Double_t               fPhiMax;                        // maximum phi
+  Double_t               fEtaMin;                        // minimum eta
+  Double_t               fEtaMax;                        // maximum eta
+  Double_t               fAreaCut;                       // cut on jet area
+  TF1                   *fScaleFunction;                 // pre-computed scale factor as a function of centrality
+  TList                 *fOutputList;                    //!output list
+  TH1F                  *fHistCentrality;                //!centrality distribution
+  TH1F                  *fHistJetPt;                     //!jet pt distribution
+  TH1F                  *fHistJetArea;                   //!jet area
+  TH2F                  *fHistRhovsCent;                 //!rho vs. centrality
+  TH2F                  *fHistDeltaRhovsCent;            //!delta rho vs. centrality
+  TH2F                  *fHistDeltaRhoScalevsCent;       //!delta rhoscaled vs. centrality
+  TH2F                  *fHistJetPtvsCent;               //!jet pt vs. centrality
+  TH2F                  *fHistJetAreavsCent;             //!jet area vs. centrality
+  TH2F                  *fHistNjetvsCent;                //!no. of jets vs. centrality
+  TH2F                  *fHistRhovsNtrack;               //!rho vs. no. of tracks
+  TH2F                  *fHistDeltaRhovsNtrack;          //!delta rho vs. no. of tracks
+  TH2F                  *fHistDeltaRhoScalevsNtrack;     //!delta rho scaled vs. no. of tracks
+  TH2F                  *fHistJetPtvsNtrack;             //!jet pt vs. no. of tracks
+  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
   
-  ClassDef(AliAnalysisTaskRho, 1); // Rho task
+  ClassDef(AliAnalysisTaskRho, 3); // Rho task
 };
-
 #endif