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   fSigJetsName("SJets")
  28 {
  29   // Constructor.
  30 }
  31
  32 //________________________________________________________________________
  33 AliAnalysisTaskRhoSparse::AliAnalysisTaskRhoSparse(const char *name, Bool_t histo) :
  34   AliAnalysisTaskRhoBase(name, histo),
  35   fHistOccCorrvsCent(0),
  36   fNExclLeadJets(0),
  37   fRhoCMS(0),
  38   fSigJetsName("SJets")
  39 {
  40   // Constructor.
  41 }
  42
  43 //________________________________________________________________________
  44 void AliAnalysisTaskRhoSparse::UserCreateOutputObjects()
  45 {
  46   if (!fCreateHisto)
  47     return;
  48
  49   AliAnalysisTaskRhoBase::UserCreateOutputObjects();
  50
  51   fHistOccCorrvsCent = new TH2F("OccCorrvsCent", "OccCorrvsCent", 101, -1, 100, 2000, 0 , 2);
  52   fOutput->Add(fHistOccCorrvsCent);
  53 }
  54
  55 //________________________________________________________________________
  56 Bool_t AliAnalysisTaskRhoSparse::IsJetOverlapping(AliEmcalJet* jet1, AliEmcalJet* jet2)
  57 {
  58   for (Int_t i = 0; i < jet1->GetNumberOfTracks(); ++i)
  59   {
  60     Int_t jet1Track = jet1->TrackAt(i);
  61     for (Int_t j = 0; j < jet2->GetNumberOfTracks(); ++j)
  62     {
  63       Int_t jet2Track = jet2->TrackAt(j);
  64       if (jet1Track == jet2Track)
  65         return kTRUE;
  66     }
  67   }
  68   return kFALSE;
  69 }
  70
  71 //________________________________________________________________________
  72 Bool_t AliAnalysisTaskRhoSparse::IsJetSignal(AliEmcalJet* jet)
  73 {
  74   if(jet->Pt()>5){
  75       return kTRUE;
  76   }else{
  77     return kFALSE;
  78   }
  79 }
  80
  81
  82 //________________________________________________________________________
  83 Bool_t AliAnalysisTaskRhoSparse::Run()
  84 {
  85   // Run the analysis.
  86
  87   fRho->SetVal(0);
  88   if (fRhoScaled)
  89     fRhoScaled->SetVal(0);
  90
  91   if (!fJets)
  92     return kFALSE;
  93
  94   const Int_t Njets   = fJets->GetEntries();
  95
  96   TClonesArray *sigjets = 0;
  97   sigjets= dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fSigJetsName));
  98
  99   Int_t NjetsSig = 0;
 100   if(sigjets) NjetsSig = sigjets->GetEntries();
 101
 102
 103   Int_t maxJetIds[]   = {-1, -1};
 104   Float_t maxJetPts[] = { 0,  0};
 105
 106   if (fNExclLeadJets > 0) {
 107     for (Int_t ij = 0; ij < Njets; ++ij) {
 108       AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
 109       if (!jet) {
 110         AliError(Form("%s: Could not receive jet %d", GetName(), ij));
 111         continue;
 112       }
 113
 114       if (!AcceptJet(jet))
 115         continue;
 116
 117       if (jet->Pt() > maxJetPts[0]) {
 118         maxJetPts[1] = maxJetPts[0];
 119         maxJetIds[1] = maxJetIds[0];
 120         maxJetPts[0] = jet->Pt();
 121         maxJetIds[0] = ij;
 122       } else if (jet->Pt() > maxJetPts[1]) {
 123         maxJetPts[1] = jet->Pt();
 124         maxJetIds[1] = ij;
 125       }
 126     }
 127     if (fNExclLeadJets < 2) {
 128       maxJetIds[1] = -1;
 129       maxJetPts[1] = 0;
 130     }
 131   }
 132
 133   static Double_t rhovec[999];
 134   Int_t NjetAcc = 0;
 135   Double_t TotaljetArea=0;
 136   Double_t TotaljetAreaPhys=0;
 137
 138   // push all jets within selected acceptance into stack
 139   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
 140
 141     // exlcuding lead jets
 142     if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
 143       continue;
 144
 145     AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
 146     if (!jet) {
 147       AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
 148       continue;
 149     }
 150
 151     TotaljetArea+=jet->Area();
 152
 153     if(jet->Pt()>0.1){
 154       TotaljetAreaPhys+=jet->Area();
 155     }
 156
 157     if (!AcceptJet(jet))
 158       continue;
 159
 160
 161
 162    // Search for overlap with signal jets
 163     Bool_t isOverlapping = kFALSE;
 164     if(sigjets){
 165       for(Int_t j=0;j<NjetsSig;j++)
 166         {
 167           AliEmcalJet* signalJet = static_cast<AliEmcalJet*>(sigjets->At(j));
 168           if(!AcceptJet(signalJet))
 169             continue;
 170           if(!IsJetSignal(signalJet))
 171             continue;
 172
 173           if(IsJetOverlapping(signalJet, jet))
 174             {
 175               isOverlapping = kTRUE;
 176               break;
 177             }
 178         }
 179     }
 180
 181     if(isOverlapping)
 182       continue;
 183
 184     if(jet->Pt()>0.1){
 185       rhovec[NjetAcc] = jet->Pt() / jet->Area();
 186       ++NjetAcc;
 187     }
 188
 189
 190   }
 191
 192   Double_t OccCorr=0.0;
 193   if(TotaljetArea>0) OccCorr=TotaljetAreaPhys/TotaljetArea;
 194
 195   fHistOccCorrvsCent->Fill(fCent, OccCorr);
 196
 197
 198   if (NjetAcc > 0) {
 199     //find median value
 200     Double_t rho = TMath::Median(NjetAcc, rhovec);
 201
 202
 203     if(fRhoCMS){
 204       rho = rho * OccCorr;
 205     }
 206
 207     fRho->SetVal(rho);
 208
 209     if (fRhoScaled) {
 210       Double_t rhoScaled = rho * GetScaleFactor(fCent);
 211       fRhoScaled->SetVal(rhoScaled);
 212     }
 213   }
 214
 215   return kTRUE;
 216 }