add tag status cut
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetContainer.cxx
1 // $Id$
2 //
3 // Container with name, TClonesArray and cuts for jets
4 //
5 // Author: M. Verweij
6
7 #include <TClonesArray.h>
8
9 #include "AliEmcalJet.h"
10 #include "AliVEvent.h"
11 #include "AliLog.h"
12 #include "AliEMCALGeometry.h"
13 #include "AliParticleContainer.h"
14 #include "AliClusterContainer.h"
15 #include "AliLocalRhoParameter.h"
16
17 #include "AliJetContainer.h"
18
19 ClassImp(AliJetContainer)
20
21 //________________________________________________________________________
22 AliJetContainer::AliJetContainer():
23   AliEmcalContainer("AliJetContainer"),
24   fJetAcceptanceType(kUser),
25   fJetRadius(0),
26   fRhoName(),
27   fLocalRhoName(),
28   fFlavourSelection(0),
29   fPtBiasJetTrack(0),
30   fPtBiasJetClus(0),
31   fJetPtCut(1),
32   fJetAreaCut(-1),
33   fAreaEmcCut(0),
34   fJetMinEta(-0.9),
35   fJetMaxEta(0.9),
36   fJetMinPhi(-10),
37   fJetMaxPhi(10),
38   fMaxClusterPt(1000),
39   fMaxTrackPt(100),
40   fZLeadingEmcCut(10.),
41   fZLeadingChCut(10.),
42   fNEFMinCut(-10.),
43   fNEFMaxCut(10.),
44   fLeadingHadronType(0),
45   fNLeadingJets(1),
46   fJetBitMap(0),
47   fJetTrigger(0),
48   fTagStatus(-1),
49   fParticleContainer(0),
50   fClusterContainer(0),
51   fRho(0),
52   fLocalRho(0),
53   fGeom(0),
54   fRunNumber(0)
55 {
56   // Default constructor.
57
58   fClassName = "AliEmcalJet";
59 }
60
61 //________________________________________________________________________
62 AliJetContainer::AliJetContainer(const char *name):
63   AliEmcalContainer(name),
64   fJetAcceptanceType(kUser),
65   fJetRadius(0),
66   fRhoName(),
67   fLocalRhoName(),
68   fFlavourSelection(0),
69   fPtBiasJetTrack(0),
70   fPtBiasJetClus(0),
71   fJetPtCut(1),
72   fJetAreaCut(-1),
73   fAreaEmcCut(0),
74   fJetMinEta(-0.9),
75   fJetMaxEta(0.9),
76   fJetMinPhi(-10),
77   fJetMaxPhi(10),
78   fMaxClusterPt(1000),
79   fMaxTrackPt(100),
80   fZLeadingEmcCut(10.),
81   fZLeadingChCut(10.),
82   fNEFMinCut(-10.),
83   fNEFMaxCut(10.),
84   fLeadingHadronType(0),
85   fNLeadingJets(1),
86   fJetBitMap(0),
87   fJetTrigger(0),
88   fTagStatus(-1),
89   fParticleContainer(0),
90   fClusterContainer(0),
91   fRho(0),
92   fLocalRho(0),
93   fGeom(0),
94   fRunNumber(0)
95 {
96   // Standard constructor.
97
98   fClassName = "AliEmcalJet";
99 }
100
101 //________________________________________________________________________
102 void AliJetContainer::SetArray(AliVEvent *event) 
103 {
104   // Set jet array
105
106   AliEmcalContainer::SetArray(event);
107
108   if(fJetAcceptanceType==kTPC) {
109     AliDebug(2,Form("%s: set TPC acceptance cuts",GetName()));
110     SetJetEtaPhiTPC();
111   }
112   else if(fJetAcceptanceType==kEMCAL) {
113     AliDebug(2,Form("%s: set EMCAL acceptance cuts",GetName()));
114     SetJetEtaPhiEMCAL();
115  }
116 }
117
118
119 //________________________________________________________________________
120 void AliJetContainer::SetEMCALGeometry() {
121   fGeom = AliEMCALGeometry::GetInstance();
122   if (!fGeom) {
123     AliError(Form("%s: Can not create geometry", GetName()));
124     return;
125   }
126 }
127
128 //________________________________________________________________________
129 void AliJetContainer::LoadRho(AliVEvent *event)
130 {
131   // Load rho
132
133   if (!fRhoName.IsNull() && !fRho) {
134     fRho = dynamic_cast<AliRhoParameter*>(event->FindListObject(fRhoName));
135     if (!fRho) {
136       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
137       return;
138     }
139   }
140 }
141
142 //________________________________________________________________________
143 void AliJetContainer::LoadLocalRho(AliVEvent *event)
144 {
145   // Load local rho
146
147   if (!fLocalRhoName.IsNull() && !fLocalRho) {
148     fLocalRho = dynamic_cast<AliLocalRhoParameter*>(event->FindListObject(fLocalRhoName));
149     if (!fLocalRho) {
150       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fLocalRhoName.Data()));
151       return;
152     }
153   }
154 }
155
156 //________________________________________________________________________
157 AliEmcalJet* AliJetContainer::GetLeadingJet(const char* opt)
158 {
159   // Get the leading jet; if opt contains "rho" the sorting is according to pt-A*rho
160
161   TString option(opt);
162   option.ToLower();
163
164   Int_t tempID = fCurrentID;
165
166   AliEmcalJet *jetMax = GetNextAcceptJet(0);
167   AliEmcalJet *jet = 0;
168
169   if (option.Contains("rho")) {
170     while ((jet = GetNextAcceptJet())) {
171       if ( (jet->Pt()-jet->Area()*GetRhoVal()) > (jetMax->Pt()-jetMax->Area()*GetRhoVal()) )
172         jetMax = jet;
173     }
174   }
175   else {
176     while ((jet = GetNextAcceptJet())) {
177       if (jet->Pt() > jetMax->Pt()) jetMax = jet;
178     }
179   }
180
181   fCurrentID = tempID;
182
183   return jetMax;
184 }
185
186 //________________________________________________________________________
187 AliEmcalJet* AliJetContainer::GetJet(Int_t i) const {
188
189   //Get i^th jet in array
190
191   if(i<0 || i>fClArray->GetEntriesFast()) return 0;
192   AliEmcalJet *jet = static_cast<AliEmcalJet*>(fClArray->At(i));
193   return jet;
194
195 }
196
197 //________________________________________________________________________
198 AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) const {
199
200   //Only return jet if is accepted
201
202   AliEmcalJet *jet = GetJet(i);
203   if(!AcceptJet(jet)) return 0;
204
205   return jet;
206 }
207
208 //________________________________________________________________________
209 AliEmcalJet* AliJetContainer::GetJetWithLabel(Int_t lab) const {
210
211   //Get particle with label lab in array
212   
213   Int_t i = GetIndexFromLabel(lab);
214   return GetJet(i);
215 }
216
217 //________________________________________________________________________
218 AliEmcalJet* AliJetContainer::GetAcceptJetWithLabel(Int_t lab) const {
219
220   //Get particle with label lab in array
221   
222   Int_t i = GetIndexFromLabel(lab);
223   return GetAcceptJet(i);
224 }
225
226 //________________________________________________________________________
227 AliEmcalJet* AliJetContainer::GetNextAcceptJet(Int_t i) {
228
229   //Get next accepted jet; if i >= 0 (re)start counter from i; return 0 if no accepted jet could be found
230
231   if (i>=0) fCurrentID = i;
232
233   const Int_t njets = GetNEntries();
234   AliEmcalJet *jet = 0;
235   while (fCurrentID < njets && !jet) { 
236     jet = GetAcceptJet(fCurrentID);
237     fCurrentID++;
238   }
239
240   return jet;
241 }
242
243 //________________________________________________________________________
244 AliEmcalJet* AliJetContainer::GetNextJet(Int_t i) {
245
246   //Get next jet; if i >= 0 (re)start counter from i; return 0 if no jet could be found
247
248   if (i>=0) fCurrentID = i;
249
250   const Int_t njets = GetNEntries();
251   AliEmcalJet *jet = 0;
252   while (fCurrentID < njets && !jet) { 
253     jet = GetJet(fCurrentID);
254     fCurrentID++;
255   }
256
257   return jet;
258 }
259
260 //________________________________________________________________________
261 Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
262   AliEmcalJet *jet = GetJet(i);
263
264   return jet->Pt() - fRho->GetVal()*jet->Area();
265 }
266
267 //________________________________________________________________________
268 Double_t AliJetContainer::GetJetPtCorrLocal(Int_t i) const {
269   AliEmcalJet *jet = GetJet(i);
270
271   return jet->Pt() - fLocalRho->GetLocalVal(jet->Phi(), fJetRadius)*jet->Area();
272 }
273
274 //________________________________________________________________________
275 void AliJetContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
276 {
277   //Get momentum of the i^th jet in array
278
279   AliEmcalJet *jet = GetJet(i);
280   if(jet) jet->GetMom(mom);
281 }
282
283 //________________________________________________________________________
284 Bool_t AliJetContainer::AcceptBiasJet(AliEmcalJet *jet) const
285
286   // Accept jet with a bias.
287
288   if (fLeadingHadronType == 0) {
289     if (jet->MaxTrackPt() < fPtBiasJetTrack) return kFALSE;
290   }
291   else if (fLeadingHadronType == 1) {
292     if (jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
293   }
294   else {
295     if (jet->MaxTrackPt() < fPtBiasJetTrack && jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
296   }
297
298   return kTRUE;
299
300
301 }
302
303 //________________________________________________________________________
304 Bool_t AliJetContainer::AcceptJet(AliEmcalJet *jet) const
305 {   
306
307    // Return true if jet is accepted.
308
309    if (!jet)
310       return kFALSE;
311
312    if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap)
313       return kFALSE;
314
315    if (jet->Pt() <= fJetPtCut) 
316       return kFALSE;
317
318    if (jet->Area() <= fJetAreaCut) 
319       return kFALSE;
320
321    if (jet->AreaEmc() < fAreaEmcCut)
322       return kFALSE;
323    
324    if (fZLeadingChCut < 1 && GetZLeadingCharged(jet) > fZLeadingChCut)
325       return kFALSE;
326    
327    if (fZLeadingEmcCut < 1 && GetZLeadingEmc(jet) > fZLeadingEmcCut)
328       return kFALSE;
329
330    if (jet->NEF() < fNEFMinCut || jet->NEF() > fNEFMaxCut)
331       return kFALSE;
332    
333    if (!AcceptBiasJet(jet))
334       return kFALSE;
335    
336    if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
337       return kFALSE;
338    
339    if (fFlavourSelection != 0 && !jet->TestFlavourTag(fFlavourSelection))
340       return kFALSE;
341
342    if(fTagStatus>-1 && jet->GetTagStatus()!=fTagStatus)
343      return kFALSE;
344    
345    Double_t jetPhi = jet->Phi();
346    Double_t jetEta = jet->Eta();
347    
348    if (fJetMinPhi < 0) // if limits are given in (-pi, pi) range
349       jetPhi -= TMath::Pi() * 2;
350    
351    return (Bool_t)(jetEta > fJetMinEta && jetEta < fJetMaxEta && jetPhi > fJetMinPhi && jetPhi < fJetMaxPhi);
352 }
353
354 //________________________________________________________________________
355 Double_t AliJetContainer::GetLeadingHadronPt(AliEmcalJet *jet) const
356 {
357   if (fLeadingHadronType == 0)       // charged leading hadron
358     return jet->MaxTrackPt();
359   else if (fLeadingHadronType == 1)  // neutral leading hadron
360     return jet->MaxClusterPt();
361   else                               // charged or neutral
362     return jet->MaxPartPt();
363 }
364
365 //________________________________________________________________________
366 void AliJetContainer::GetLeadingHadronMomentum(TLorentzVector &mom, AliEmcalJet *jet) const
367 {
368   Double_t maxClusterPt = 0;
369   Double_t maxClusterEta = 0;
370   Double_t maxClusterPhi = 0;
371   
372   Double_t maxTrackPt = 0;
373   Double_t maxTrackEta = 0;
374   Double_t maxTrackPhi = 0;
375       
376   if (fClusterContainer && fClusterContainer->GetArray() && (fLeadingHadronType == 1 || fLeadingHadronType == 2)) {
377     AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
378     if (cluster) {
379       TLorentzVector nPart;
380       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
381       
382       maxClusterEta = nPart.Eta();
383       maxClusterPhi = nPart.Phi();
384       maxClusterPt = nPart.Pt();
385     }
386   }
387       
388   if (fParticleContainer && fParticleContainer->GetArray() && (fLeadingHadronType == 0 || fLeadingHadronType == 2)) {
389     AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
390     if (track) {
391       maxTrackEta = track->Eta();
392       maxTrackPhi = track->Phi();
393       maxTrackPt = track->Pt();
394     }
395   }
396       
397   if (maxTrackPt > maxClusterPt) 
398     mom.SetPtEtaPhiM(maxTrackPt,maxTrackEta,maxTrackPhi,0.139);
399   else 
400     mom.SetPtEtaPhiM(maxClusterPt,maxClusterEta,maxClusterPhi,0.139);
401 }
402
403 //________________________________________________________________________
404 Double_t AliJetContainer::GetZLeadingEmc(AliEmcalJet *jet) const
405 {
406
407   if (fClusterContainer && fClusterContainer->GetArray()) {
408     TLorentzVector mom;
409     
410     AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
411     if (cluster) {
412       cluster->GetMomentum(mom, const_cast<Double_t*>(fVertex));
413       
414       return GetZ(jet,mom);
415     }
416     else
417       return -1;
418   }
419   else
420     return -1;
421 }
422
423 //________________________________________________________________________
424 Double_t AliJetContainer::GetZLeadingCharged(AliEmcalJet *jet) const
425 {
426
427   if (fParticleContainer && fParticleContainer->GetArray() ) {
428     TLorentzVector mom;
429     
430     AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
431     if (track) {
432       mom.SetPtEtaPhiM(track->Pt(),track->Eta(),track->Phi(),0.139);
433       
434       return GetZ(jet,mom);
435     }
436     else
437       return -1;
438   }
439   else
440     return -1;
441 }
442
443 //________________________________________________________________________
444 Double_t AliJetContainer::GetZ(AliEmcalJet *jet, TLorentzVector mom) const
445 {
446
447   Double_t pJetSq = jet->Px()*jet->Px() + jet->Py()*jet->Py() + jet->Pz()*jet->Pz();
448
449   if(pJetSq>1e-6)
450     return (mom.Px()*jet->Px() + mom.Py()*jet->Py() + mom.Pz()*jet->Pz())/pJetSq;
451   else {
452     AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
453     return -1;
454   }
455
456 }
457
458 //________________________________________________________________________
459 void AliJetContainer::SetJetEtaPhiEMCAL()
460 {
461   //Set default cuts for full jets
462
463   if(!fGeom) SetEMCALGeometry();
464   if(fGeom) {
465     SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius);
466
467     if(fRunNumber>=177295 && fRunNumber<=197470) //small SM masked in 2012 and 2013
468       SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
469     else
470       SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);
471
472   }
473   else {
474     AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
475     SetJetEtaLimits(-0.7+fJetRadius,0.7-fJetRadius);
476     SetJetPhiLimits(1.405+fJetRadius,3.135-fJetRadius);
477   }
478 }
479
480 //________________________________________________________________________
481 void AliJetContainer::SetJetEtaPhiTPC()
482 {
483   //Set default cuts for charged jets
484
485   SetJetEtaLimits(-0.9+fJetRadius, 0.9-fJetRadius);
486   SetJetPhiLimits(-10, 10);
487 }
488
489 //________________________________________________________________________
490 void AliJetContainer::ResetCuts() 
491 {
492   // Reset cuts to default values
493
494   fPtBiasJetTrack = 0;
495   fPtBiasJetClus = 0;
496   fJetPtCut = 1;
497   fJetAreaCut = -1;
498   fAreaEmcCut = 0;
499   fJetMinEta = -0.9;
500   fJetMaxEta = 0.9;
501   fJetMinPhi = -10;
502   fJetMaxPhi = 10;
503   fMaxClusterPt = 1000;
504   fMaxTrackPt = 100;
505   fLeadingHadronType = 0;
506   fZLeadingEmcCut = 10.;
507   fZLeadingChCut = 10.;
508 }
509
510 //________________________________________________________________________
511 void AliJetContainer::SetClassName(const char *clname)
512 {
513   // Set the class name
514
515   TClass cls(clname);
516   if (cls.InheritsFrom("AliEmcalJet")) fClassName = clname;
517   else AliError(Form("Unable to set class name %s for a AliJetContainer, it must inherits from AliEmcalJet!",clname));
518 }
519
520 //________________________________________________________________________
521 Double_t AliJetContainer::GetFractionSharedPt(AliEmcalJet *jet1) const
522 {
523   //
524   // Get fraction of shared pT between matched full and charged jet
525   // Uses charged jet pT as baseline: fraction = \Sum_{const,full jet} pT,const,i / pT,jet,ch
526   // Only works if tracks array of both jets is the same
527   //
528
529   AliEmcalJet *jet2 = jet1->ClosestJet();
530   if(!jet2) return -1;
531
532   Double_t fraction = 0.;
533   Double_t jetPt2 = jet2->Pt();
534  
535   if(jetPt2>0) {
536     Double_t sumPt = 0.;
537     AliVParticle *vpf = 0x0;
538     Int_t iFound = 0;
539     for(Int_t icc=0; icc<jet2->GetNumberOfTracks(); icc++) {
540       Int_t idx = (Int_t)jet2->TrackAt(icc);
541       iFound = 0;
542       for(Int_t icf=0; icf<jet1->GetNumberOfTracks(); icf++) {
543         if(idx == jet1->TrackAt(icf) && iFound==0 ) {
544           iFound=1;
545           vpf = static_cast<AliVParticle*>(jet1->TrackAt(icf, fParticleContainer->GetArray()));
546           if(vpf) sumPt += vpf->Pt();
547           continue;
548         }
549       }
550     }
551     fraction = sumPt/jetPt2;
552   } else 
553     fraction = -1;
554   
555   return fraction;
556 }
557