PWGJE/EMCALJetTasks/AliAnalysisTaskRhoSparse.cxx

   1 // $Id: AliAnalysisTaskRho.cxx 58408 2012-09-03 07:00:58Z loizides $
   2 //
   3 // Calculation of rho from a collection of jets.
   4 // If scale function is given the scaled rho will be exported
   5 // with the name as "fRhoName".apppend("_Scaled").
   6 //
   7 // Authors: R.Reed, S.Aiola, M.Connors
   8
   9 #include "AliAnalysisTaskRhoSparse.h"
  10
  11 #include <TClonesArray.h>
  12 #include <TMath.h>
  13
  14 #include "AliAnalysisManager.h"
  15 #include "AliEmcalJet.h"
  16 #include "AliLog.h"
  17 #include "AliRhoParameter.h"
  18
  19 ClassImp(AliAnalysisTaskRhoSparse)
  20
  21 //________________________________________________________________________
  22 AliAnalysisTaskRhoSparse::AliAnalysisTaskRhoSparse() :
  23   AliAnalysisTaskRhoBase("AliAnalysisTaskRhoSparse"),
  24   fHistOccCorrvsCent(0),
  25   fNExclLeadJets(0),
  26   fRhoCMS(0)
  27 {
  28   // Constructor.
  29 }
  30
  31 //________________________________________________________________________
  32 AliAnalysisTaskRhoSparse::AliAnalysisTaskRhoSparse(const char *name, Bool_t histo) :
  33   AliAnalysisTaskRhoBase(name, histo),
  34   fHistOccCorrvsCent(0),
  35   fNExclLeadJets(0),
  36   fRhoCMS(0)
  37 {
  38   // Constructor.
  39 }
  40
  41 //________________________________________________________________________
  42 void AliAnalysisTaskRhoSparse::UserCreateOutputObjects()
  43 {
  44   if (!fCreateHisto)
  45     return;
  46
  47   AliAnalysisTaskRhoBase::UserCreateOutputObjects();
  48
  49   fHistOccCorrvsCent = new TH2F("OccCorrvsCent", "OccCorrvsCent", 101, -1, 100, 2000, 0 , 2);
  50   fOutput->Add(fHistOccCorrvsCent);
  51 }
  52
  53 //________________________________________________________________________
  54 Bool_t AliAnalysisTaskRhoSparse::Run()
  55 {
  56   // Run the analysis.
  57
  58   fRho->SetVal(0);
  59   if (fRhoScaled)
  60     fRhoScaled->SetVal(0);
  61
  62   if (!fJets)
  63     return kFALSE;
  64
  65   const Int_t Njets   = fJets->GetEntries();
  66
  67   Int_t maxJetIds[]   = {-1, -1};
  68   Float_t maxJetPts[] = { 0,  0};
  69
  70   if (fNExclLeadJets > 0) {
  71     for (Int_t ij = 0; ij < Njets; ++ij) {
  72       AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
  73       if (!jet) {
  74         AliError(Form("%s: Could not receive jet %d", GetName(), ij));
  75         continue;
  76       }
  77
  78       if (!AcceptJet(jet))
  79         continue;
  80
  81       if (jet->Pt() > maxJetPts[0]) {
  82         maxJetPts[1] = maxJetPts[0];
  83         maxJetIds[1] = maxJetIds[0];
  84         maxJetPts[0] = jet->Pt();
  85         maxJetIds[0] = ij;
  86       } else if (jet->Pt() > maxJetPts[1]) {
  87         maxJetPts[1] = jet->Pt();
  88         maxJetIds[1] = ij;
  89       }
  90     }
  91     if (fNExclLeadJets < 2) {
  92       maxJetIds[1] = -1;
  93       maxJetPts[1] = 0;
  94     }
  95   }
  96
  97   static Double_t rhovec[999];
  98   Int_t NjetAcc = 0;
  99   Double_t TotaljetArea=0;
 100   Double_t TotaljetAreaPhys=0;
 101
 102   // push all jets within selected acceptance into stack
 103   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
 104
 105     // exlcuding lead jets
 106     if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
 107       continue;
 108
 109     AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
 110     if (!jet) {
 111       AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
 112       continue;
 113     }
 114
 115     //cout << "jetpt: " << jet->Pt() << " jetArea: " << jet->Area() <<endl;
 116
 117     if (!AcceptJet(jet))
 118       continue;
 119
 120
 121     if(jet->Pt()>0.01) rhovec[NjetAcc] = jet->Pt() / jet->Area();
 122     //cout << "ACCEPTED: jetpt: " << jet->Pt() << " jetArea: " << jet->Area() << " jetRho: " << rhovec[NjetAcc] <<endl;
 123     if(jet->Pt()>0.01) TotaljetAreaPhys+=jet->Area();
 124     TotaljetArea+=jet->Area();
 125     ++NjetAcc;
 126   }
 127
 128   const Double_t TpcMaxPhi = TMath::Pi()*2.;
 129
 130   const Double_t TpcArea = TpcMaxPhi * 2.*(0.7);
 131   Double_t OccCorr=0.0;
 132   //cout << "Area Physical: " << TotaljetAreaPhys << " total: " << TotaljetArea <<endl;
 133   //if(TotaljetArea>0) OccCorr=TotaljetAreaPhys/TotaljetArea;
 134   if(TpcArea>0) OccCorr=TotaljetAreaPhys/TpcArea;
 135
 136   fHistOccCorrvsCent->Fill(fCent, OccCorr);
 137
 138
 139   if (NjetAcc > 0) {
 140     //find median value
 141     Double_t rho = TMath::Median(NjetAcc, rhovec);
 142
 143     if(fRhoCMS){
 144       rho = rho * OccCorr;
 145     }
 146
 147     fRho->SetVal(rho);
 148
 149
 150     if (fRhoScaled) {
 151       Double_t rhoScaled = rho * GetScaleFactor(fCent);
 152       fRhoScaled->SetVal(rhoScaled);
 153     }
 154   }
 155
 156   return kTRUE;
 157 }