jet framework from Marta/Salvatore
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetContainer.cxx
1 //
2 // Jet Container
3 //
4 // Author: M. Verweij
5
6 #include <TROOT.h>
7 #include <TSystem.h>
8 #include <TInterpreter.h>
9
10 #include <TChain.h>
11 #include <TClonesArray.h>
12 #include <TList.h>
13 #include <TObject.h>
14 #include "AliEmcalJet.h"
15 #include "AliVEvent.h"
16 #include "AliLog.h"
17 #include "AliRhoParameter.h"
18 #include "AliEMCALGeometry.h"
19
20 #include "AliJetContainer.h"
21
22 ClassImp(AliJetContainer)
23
24 //________________________________________________________________________
25 AliJetContainer::AliJetContainer():
26   AliEmcalContainer("AliJetContainer"),
27   fJetAcceptanceType(kUser),
28   fJetRadius(0),
29   fRhoName(),
30   fPtBiasJetTrack(0),
31   fPtBiasJetClus(0),
32   fJetPtCut(1),
33   fJetAreaCut(-1),
34   fAreaEmcCut(0),
35   fJetMinEta(-0.9),
36   fJetMaxEta(0.9),
37   fJetMinPhi(-10),
38   fJetMaxPhi(10),
39   fMaxClusterPt(1000),
40   fMaxTrackPt(100),
41   fLeadingHadronType(0),
42   fNLeadingJets(1),
43   fJetBitMap(0),
44   fRho(0),
45   fGeom(0)
46 {
47   // Default constructor.
48
49 }
50
51 //________________________________________________________________________
52 AliJetContainer::AliJetContainer(const char *name):
53   AliEmcalContainer(name),
54   fJetAcceptanceType(kUser),
55   fJetRadius(0),
56   fRhoName(),
57   fPtBiasJetTrack(0),
58   fPtBiasJetClus(0),
59   fJetPtCut(1),
60   fJetAreaCut(-1),
61   fAreaEmcCut(0),
62   fJetMinEta(-0.9),
63   fJetMaxEta(0.9),
64   fJetMinPhi(-10),
65   fJetMaxPhi(10),
66   fMaxClusterPt(1000),
67   fMaxTrackPt(100),
68   fLeadingHadronType(0),
69   fNLeadingJets(1),
70   fJetBitMap(0),
71   fRho(0),
72   fGeom(0)
73 {
74   // Standard constructor.
75
76 }
77
78 //________________________________________________________________________
79 AliJetContainer::~AliJetContainer()
80 {
81   // Destructor.
82 }
83
84 //________________________________________________________________________
85 void AliJetContainer::SetJetArray(AliVEvent *event) 
86 {
87   // Set jet array
88
89   //  SetArray(event, fClArrayName.GetName(), "AliEmcalJet");
90   SetArray(event, "AliEmcalJet");
91
92   if(fJetAcceptanceType==kTPC)
93     SetJetEtaPhiTPC();
94   else if(fJetAcceptanceType==kEMCAL)
95     SetJetEtaPhiEMCAL();
96 }
97
98 //________________________________________________________________________
99 void AliJetContainer::SetEMCALGeometry() {
100   fGeom = AliEMCALGeometry::GetInstance();
101   if (!fGeom) {
102     AliError(Form("%s: Can not create geometry", GetName()));
103     return;
104   }
105 }
106
107 //________________________________________________________________________
108 void AliJetContainer::LoadRho(AliVEvent *event)
109 {
110   //Load rho
111
112   if (!fRhoName.IsNull() && !fRho) {
113     fRho = dynamic_cast<AliRhoParameter*>(event->FindListObject(fRhoName));
114     if (!fRho) {
115       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
116       return;
117     }
118   }
119 }
120
121 //________________________________________________________________________
122 AliEmcalJet* AliJetContainer::GetJet(Int_t i) const {
123
124   //Get i^th jet in array
125
126   if(i<0 || i>fClArray->GetEntriesFast()) return 0;
127   AliEmcalJet *jet = static_cast<AliEmcalJet*>(fClArray->At(i));
128   return jet;
129
130 }
131
132
133 //________________________________________________________________________
134 Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
135   AliEmcalJet *jet = GetJet(i);
136
137   return jet->Pt() - fRho->GetVal()*jet->Area();
138 }
139
140 //________________________________________________________________________
141 AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) const {
142
143   //Only return jet if is accepted
144
145   AliEmcalJet *jet = GetJet(i);
146   if(!AcceptJet(jet)) return 0;
147
148   return jet;
149 }
150
151 //________________________________________________________________________
152 Bool_t AliJetContainer::AcceptBiasJet(AliEmcalJet *jet) const
153
154   // Accept jet with a bias.
155
156   if (fLeadingHadronType == 0) {
157     if (jet->MaxTrackPt() < fPtBiasJetTrack) return kFALSE;
158   }
159   else if (fLeadingHadronType == 1) {
160     if (jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
161   }
162   else {
163     if (jet->MaxTrackPt() < fPtBiasJetTrack && jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
164   }
165
166   return kTRUE;
167
168
169 }
170
171 //________________________________________________________________________
172 Bool_t AliJetContainer::AcceptJet(AliEmcalJet *jet) const
173 {   
174   // Return true if jet is accepted.
175   if (!jet)
176     return kFALSE;
177
178   if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap)
179     return kFALSE;
180
181   if (jet->Pt() <= fJetPtCut) 
182     return kFALSE;
183  
184   if (jet->Area() <= fJetAreaCut) 
185     return kFALSE;
186
187   if (jet->AreaEmc()<fAreaEmcCut)
188     return kFALSE;
189
190   if (!AcceptBiasJet(jet))
191     return kFALSE;
192   
193   if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
194     return kFALSE;
195
196   Double_t jetPhi = jet->Phi();
197   Double_t jetEta = jet->Eta();
198   
199   if (fJetMinPhi < 0) // if limits are given in (-pi, pi) range
200     jetPhi -= TMath::Pi() * 2;
201
202   return (Bool_t)(jetEta > fJetMinEta && jetEta < fJetMaxEta && jetPhi > fJetMinPhi && jetPhi < fJetMaxPhi);
203 }
204
205 //________________________________________________________________________
206 Double_t AliJetContainer::GetLeadingHadronPt(AliEmcalJet *jet) const
207 {
208   if (fLeadingHadronType == 0)       // charged leading hadron
209     return jet->MaxTrackPt();
210   else if (fLeadingHadronType == 1)  // neutral leading hadron
211     return jet->MaxClusterPt();
212   else                               // charged or neutral
213     return jet->MaxPartPt();
214 }
215
216
217 //________________________________________________________________________
218 void AliJetContainer::SetJetEtaPhiEMCAL()
219 {
220   //Set default cuts for full jets
221
222   if(!fGeom) SetEMCALGeometry();
223   if(fGeom) {
224     SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius);
225     SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);
226   }
227   else {
228     AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings");
229     SetJetEtaLimits(-0.7+fJetRadius,0.7-fJetRadius);
230     SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
231   }
232
233 }
234
235 //________________________________________________________________________
236 void AliJetContainer::SetJetEtaPhiTPC()
237 {
238   //Set default cuts for full jets
239
240   SetJetEtaLimits(-0.5, 0.5);
241   SetJetPhiLimits(-10, 10);
242
243 }
244
245 //________________________________________________________________________
246 void AliJetContainer::ResetCuts() 
247 {
248   // Reset cuts to default values
249
250   fPtBiasJetTrack = 0;
251   fPtBiasJetClus = 0;
252   fJetPtCut = 1;
253   fJetAreaCut = -1;
254   fAreaEmcCut = 0;
255   fJetMinEta = -0.9;
256   fJetMaxEta = 0.9;
257   fJetMinPhi = -10;
258   fJetMaxPhi = 10;
259   fMaxClusterPt = 1000;
260   fMaxTrackPt = 100;
261   fLeadingHadronType = 0;
262
263 }