setter to assume pion mass for clusters
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliAnalysisTaskRhoSparse.cxx
1 // $Id$
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   if (fCreateHisto)
188     fHistOccCorrvsCent->Fill(fCent, OccCorr);
189
190   if (NjetAcc > 0) {
191     //find median value
192     Double_t rho = TMath::Median(NjetAcc, rhovec);
193
194     if(fRhoCMS){
195       rho = rho * OccCorr;
196     }
197
198     fOutRho->SetVal(rho);
199
200     if (fOutRhoScaled) {
201       Double_t rhoScaled = rho * GetScaleFactor(fCent);
202       fOutRhoScaled->SetVal(rhoScaled);
203     }
204   }
205
206   return kTRUE;
207