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