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