4bb637582f68e67838afa16099f10d740a5c5d90
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliEmcalJetTask.cxx
1 // $Id$
2 //
3 // Emcal jet finder task.
4 //
5 // Authors: C.Loizides, S.Aiola
6
7 #include <vector>
8 #include "AliEmcalJetTask.h"
9
10 #include <TChain.h>
11 #include <TClonesArray.h>
12 #include <TList.h>
13 #include <TLorentzVector.h>
14 #include <TMath.h>
15 #include <TRandom3.h>
16
17 #include "AliAnalysisManager.h"
18 #include "AliCentrality.h"
19 #include "AliEMCALGeometry.h"
20 #include "AliESDEvent.h"
21 #include "AliEmcalJet.h"
22 #include "AliEmcalParticle.h"
23 #include "AliFJWrapper.h"
24 #include "AliMCEvent.h"
25 #include "AliVCluster.h"
26 #include "AliVEvent.h"
27 #include "AliVParticle.h"
28 using std::cout;
29 using std::endl;
30 using std::cerr;
31
32 ClassImp(AliEmcalJetTask)
33
34 //________________________________________________________________________
35 AliEmcalJetTask::AliEmcalJetTask() : 
36   AliAnalysisTaskSE("AliEmcalJetTask"),
37   fTracksName("Tracks"),
38   fCaloName("CaloClusters"),
39   fJetsName("Jets"),
40   fJetType(kNone),
41   fConstSel(0),
42   fMCConstSel(0),
43   fMarkConst(kFALSE),
44   fRadius(0.4),
45   fMinJetTrackPt(0.15),
46   fMinJetClusPt(0.15),
47   fPhiMin(0),
48   fPhiMax(TMath::TwoPi()),
49   fEtaMin(-0.9),
50   fEtaMax(+0.9),
51   fMinJetArea(0.005),
52   fMinJetPt(1.0),
53   fJetPhiMin(-10),
54   fJetPhiMax(+10),
55   fJetEtaMin(-1),
56   fJetEtaMax(+1),
57   fGhostArea(0.005),
58   fMinMCLabel(0),
59   fRecombScheme(-1),
60   fTrackEfficiency(1.),
61   fIsInit(0),
62   fIsPSelSet(0),
63   fIsMcPart(0),
64   fIsEmcPart(0),
65   fLegacyMode(kFALSE),
66   fJets(0),
67   fEvent(0),
68   fTracks(0),
69   fClus(0)
70 {
71   // Default constructor.
72 }
73
74 //________________________________________________________________________
75 AliEmcalJetTask::AliEmcalJetTask(const char *name) : 
76   AliAnalysisTaskSE(name),
77   fTracksName("Tracks"),
78   fCaloName("CaloClusters"),
79   fJetsName("Jets"),
80   fJetType(kAKT|kFullJet|kRX1Jet),
81   fConstSel(0),
82   fMCConstSel(0),
83   fMarkConst(kFALSE),
84   fRadius(0.4),
85   fMinJetTrackPt(0.15),
86   fMinJetClusPt(0.15),
87   fPhiMin(0),
88   fPhiMax(TMath::TwoPi()),
89   fEtaMin(-0.9),
90   fEtaMax(+0.9),
91   fMinJetArea(0.001),
92   fMinJetPt(1.0),
93   fJetPhiMin(-10),
94   fJetPhiMax(+10),
95   fJetEtaMin(-1),
96   fJetEtaMax(+1),
97   fGhostArea(0.005),
98   fMinMCLabel(0),
99   fRecombScheme(-1),
100   fTrackEfficiency(1.),
101   fIsInit(0),
102   fIsPSelSet(0),
103   fIsMcPart(0),
104   fIsEmcPart(0),
105   fLegacyMode(false),
106   fJets(0),
107   fEvent(0),
108   fTracks(0),
109   fClus(0)
110 {
111   // Standard constructor.
112
113   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
114 }
115
116 //________________________________________________________________________
117 AliEmcalJetTask::~AliEmcalJetTask()
118 {
119   // Destructor
120 }
121
122 //________________________________________________________________________
123 void AliEmcalJetTask::UserCreateOutputObjects()
124 {
125   // Create user objects.
126
127   fJets = new TClonesArray("AliEmcalJet");
128   fJets->SetName(fJetsName);
129 }
130
131 //________________________________________________________________________
132 void AliEmcalJetTask::UserExec(Option_t *) 
133 {
134   // Main loop, called for each event.
135   if (!fIsInit) {
136     if (!DoInit())
137       return;
138     fIsInit = kTRUE;
139   }
140
141   // clear the jet array (normally a null operation)
142   fJets->Delete();
143
144   FindJets();
145 }
146
147 //________________________________________________________________________
148 void AliEmcalJetTask::Terminate(Option_t *) 
149 {
150   // Called once at the end of the analysis.
151 }
152
153 //________________________________________________________________________
154 void AliEmcalJetTask::FindJets()
155 {
156   // Find jets.
157   if (!fTracks && !fClus){
158     cout << "WARNING NO TRACKS OR CLUSTERS:"  <<endl;
159     return;
160   }
161
162   TString name("kt");
163   fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
164   if ((fJetType & kAKT) != 0) {
165     name  = "antikt";
166     jalgo = fastjet::antikt_algorithm;
167     AliDebug(1,"Using AKT algorithm");
168   }
169   else {
170     AliDebug(1,"Using KT algorithm");
171   }
172
173   if ((fJetType & kR020Jet) != 0)
174     fRadius = 0.2;
175   else if ((fJetType & kR030Jet) != 0)
176     fRadius = 0.3; 
177   else if ((fJetType & kR040Jet) != 0)
178     fRadius = 0.4;
179
180   // setup fj wrapper
181   AliFJWrapper fjw(name, name);
182   fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
183   fjw.SetGhostArea(fGhostArea);
184   fjw.SetR(fRadius);
185   fjw.SetAlgorithm(jalgo);
186   if(fRecombScheme>0)
187     fjw.SetRecombScheme(static_cast<fastjet::RecombinationScheme>(fRecombScheme)); 
188   fjw.SetMaxRap(fEtaMax);
189   fjw.Clear();
190
191   // get primary vertex
192   Double_t vertex[3] = {0, 0, 0};
193   fEvent->GetPrimaryVertex()->GetXYZ(vertex);
194
195   AliDebug(2,Form("Jet type = %d", fJetType));
196
197   if ((fIsMcPart || ((fJetType & kFullJet) != 0) || ((fJetType & kChargedJet) != 0)) && fTracks) {
198     const Int_t Ntracks = fTracks->GetEntries();
199     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
200       AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(iTracks));
201       if (!t)
202         continue;
203       if (fIsMcPart) {
204         if (((fJetType & kChargedJet) != 0) && (t->Charge() == 0)) {
205           AliDebug(2,Form("Skipping track %d because it is neutral.", iTracks));
206           continue;
207         }
208         if (((fJetType & kNeutralJet) != 0) && (t->Charge() != 0)) {
209           AliDebug(2,Form("Skipping track %d because it is charged.", iTracks));
210           continue;
211         }
212       }
213       if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) {
214         if (t->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
215           AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask)));
216           continue;
217         }
218         else {
219           AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask)));
220         }
221       }
222       else {
223         if (t->TestBits(fConstSel) != (Int_t)fConstSel) {
224           AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));
225           continue;
226         }
227         else {
228           AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));           
229         }
230       }     
231       if (t->Pt() < fMinJetTrackPt) 
232         continue;
233       Double_t eta = t->Eta();
234       Double_t phi = t->Phi();
235       if ((eta<fEtaMin) || (eta>fEtaMax) ||
236           (phi<fPhiMin) || (phi>fPhiMax))
237         continue;
238
239       // artificial inefficiency
240       if (fTrackEfficiency < 1.) {
241         Double_t rnd = gRandom->Rndm();
242         if (fTrackEfficiency < rnd) {
243           AliDebug(2,Form("Track %d rejected due to artificial tracking inefficiency", iTracks));
244           continue;
245         }
246       }
247
248       // offset of 100 for consistency with cluster ids
249       AliDebug(2,Form("Track %d accepted (label = %d, pt = %f)", iTracks, t->GetLabel(), t->Pt()));
250       fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);  
251     }
252   }
253
254   if ((((fJetType & kFullJet) != 0) || ((fJetType & kNeutralJet) != 0)) && fClus) {
255     const Int_t Nclus = fClus->GetEntries();
256     for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
257       AliVCluster *c = 0;
258       Double_t cEta=0,cPhi=0,cPt=0;
259       Double_t cPx=0,cPy=0,cPz=0;
260       if (fIsEmcPart) {
261         AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(iClus));
262         if (!ep)
263           continue;
264
265         c = ep->GetCluster();
266         if (!c)
267           continue;
268
269         if (c->GetLabel() > fMinMCLabel) {
270           if (ep->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
271             AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
272             continue;
273           }
274           else {
275             AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
276           }
277         }
278         else {
279           if (ep->TestBits(fConstSel) != (Int_t)fConstSel) {
280             AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));
281             continue;
282           }
283           else {
284             AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));        
285           }
286         }
287
288         cEta = ep->Eta();
289         cPhi = ep->Phi();
290         cPt  = ep->Pt();
291         cPx  = ep->Px();
292         cPy  = ep->Py();
293         cPz  = ep->Pz();
294       } else {
295         c = static_cast<AliVCluster*>(fClus->At(iClus));
296         if (!c)
297           continue;
298
299         if (c->GetLabel() > fMinMCLabel) {
300           if (c->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
301             AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
302             continue;
303           }
304           else {
305             AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
306           }
307         }
308         else {
309           if (c->TestBits(fConstSel) != (Int_t)fConstSel) {
310             AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));
311             continue;
312           }
313           else {
314             AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));         
315           }
316         }
317
318         TLorentzVector nP;
319         c->GetMomentum(nP, vertex);
320         cEta = nP.Eta();
321         cPhi = nP.Phi();
322         cPt  = nP.Pt();
323         cPx  = nP.Px();
324         cPy  = nP.Py();
325         cPz  = nP.Pz();
326       }
327       if (!c->IsEMCAL())
328         continue;
329       if (cPt < fMinJetClusPt) 
330         continue;
331       if ((cEta<fEtaMin) || (cEta>fEtaMax) ||
332           (cPhi<fPhiMin) || (cPhi>fPhiMax))
333         continue;
334       // offset of 100 to skip ghost particles uid = -1
335       AliDebug(2,Form("Cluster %d accepted (label = %d)", iClus, c->GetLabel()));
336       fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz), -iClus - 100);
337     }
338   }
339
340   // setting legacy mode
341   if (fLegacyMode) { 
342     fjw.SetLegacyMode(kTRUE);
343   }
344
345   // run jet finder
346   fjw.Run();
347
348   // get geometry
349   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
350   if (!geom) {
351     AliFatal(Form("%s: Can not create geometry", GetName()));
352     return;
353   }
354
355   // loop over fastjet jets
356   std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
357   // sort jets according to jet pt
358   static Int_t indexes[9999] = {-1};
359   GetSortedArray(indexes, jets_incl);
360
361   AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
362   for (UInt_t ijet=0, jetCount=0; ijet<jets_incl.size(); ++ijet) {
363     Int_t ij = indexes[ijet];
364     AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fjw.GetJetArea(ij)));
365
366     if (jets_incl[ij].perp()<fMinJetPt) 
367       continue;
368     if (fjw.GetJetArea(ij)<fMinJetArea)
369       continue;
370     if ((jets_incl[ij].eta()<fJetEtaMin) || (jets_incl[ij].eta()>fJetEtaMax) ||
371         (jets_incl[ij].phi()<fJetPhiMin) || (jets_incl[ij].phi()>fJetPhiMax))
372       continue;
373
374     AliEmcalJet *jet = new ((*fJets)[jetCount]) 
375       AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
376
377     // loop over constituents
378     std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
379     jet->SetNumberOfTracks(constituents.size());
380     jet->SetNumberOfClusters(constituents.size());
381
382     Int_t nt            = 0;
383     Int_t nc            = 0;
384     Double_t neutralE   = 0;
385     Double_t maxCh      = 0;
386     Double_t maxNe      = 0;
387     Int_t gall          = 0;
388     Int_t gemc          = 0;
389     Int_t cemc          = 0;
390     Int_t ncharged      = 0;
391     Int_t nneutral      = 0;
392     Double_t mcpt       = 0;
393     Double_t emcpt      = 0;
394
395     for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
396       Int_t uid = constituents[ic].user_index();
397       AliDebug(3,Form("Processing constituent %d", uid));
398       if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
399         ++gall;
400         Double_t gphi = constituents[ic].phi();
401         if (gphi<0) 
402           gphi += TMath::TwoPi();
403         gphi *= TMath::RadToDeg();
404         Double_t geta = constituents[ic].eta();
405         if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
406             (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
407           ++gemc; 
408       } else if ((uid > 0) && fTracks) { // track constituent
409         Int_t tid = uid - 100;
410         AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
411         if (!t) {
412           AliError(Form("Could not find track %d",tid));
413           continue;
414         }
415         if (jetCount < fMarkConst) {
416           AliDebug(2,Form("Marking track %d with bit map %d", tid, fJetType));
417           t->SetBit(fJetType);
418         }
419         Double_t cEta = t->Eta();
420         Double_t cPhi = t->Phi();
421         Double_t cPt  = t->Pt();
422         Double_t cP   = t->P();
423         if (t->Charge() == 0) {
424           neutralE += cP;
425           ++nneutral;
426           if (cPt > maxNe)
427             maxNe = cPt;
428         } else {
429           ++ncharged;
430           if (cPt > maxCh)
431             maxCh = cPt;
432         }
433         if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) // check if MC particle
434           mcpt += cPt;
435
436         if (cPhi<0) 
437           cPhi += TMath::TwoPi();
438         cPhi *= TMath::RadToDeg();
439         if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
440             (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
441           emcpt += cPt;
442           ++cemc;
443         }
444
445         jet->AddTrackAt(tid, nt);
446         ++nt;
447       } else if (fClus) { // cluster constituent
448         Int_t cid = -uid - 100;
449         AliVCluster *c = 0;
450         Double_t cEta=0,cPhi=0,cPt=0,cP=0;
451         if (fIsEmcPart) {
452           AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
453           if (!ep)
454             continue;
455           c = ep->GetCluster();
456           if (!c)
457             continue;
458           if (jetCount < fMarkConst)
459             ep->SetBit(fJetType);
460           cEta = ep->Eta();
461           cPhi = ep->Phi();
462           cPt  = ep->Pt();
463           cP   = ep->P();
464         } else {
465           c = static_cast<AliVCluster*>(fClus->At(cid));
466           if (!c)
467             continue;
468           if (jetCount < fMarkConst)
469             c->SetBit(fJetType);
470           TLorentzVector nP;
471           c->GetMomentum(nP, vertex);
472           cEta = nP.Eta();
473           cPhi = nP.Phi();
474           cPt  = nP.Pt();
475           cP   = nP.P();
476         }
477
478         neutralE += cP;
479         if (cPt > maxNe)
480           maxNe = cPt;
481
482         if (c->GetLabel() > fMinMCLabel) // MC particle
483           mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
484
485         if (cPhi<0) 
486           cPhi += TMath::TwoPi();
487         cPhi *= TMath::RadToDeg();
488         if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
489             (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
490           emcpt += cPt;
491           ++cemc;
492         }
493
494         jet->AddClusterAt(cid, nc);
495         ++nc;
496         ++nneutral;
497       } else {
498         AliError(Form("%s: No logical way to end up here.", GetName()));
499         continue;
500       }
501     } /* end of constituent loop */
502
503     jet->SetNumberOfTracks(nt);
504     jet->SetNumberOfClusters(nc);
505     jet->SortConstituents();
506     jet->SetMaxNeutralPt(maxNe);
507     jet->SetMaxChargedPt(maxCh);
508     jet->SetNEF(neutralE / jet->E());
509     fastjet::PseudoJet area(fjw.GetJetAreaVector(ij));
510     jet->SetArea(area.perp());
511     jet->SetAreaEta(area.eta());
512     jet->SetAreaPhi(area.phi());
513     jet->SetNumberOfCharged(ncharged);
514     jet->SetNumberOfNeutrals(nneutral);
515     jet->SetMCPt(mcpt);
516     jet->SetNEmc(cemc);
517     jet->SetPtEmc(emcpt);
518
519     if (gall > 0)
520       jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
521     else 
522       jet->SetAreaEmc(-1);
523     if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) && 
524         (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
525         (jet->Eta() > geom->GetArm1EtaMin()) && 
526         (jet->Eta() < geom->GetArm1EtaMax()))
527       jet->SetAxisInEmcal(kTRUE);
528
529     AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), (Int_t)constituents.size()));
530     jetCount++;
531   }
532   //fJets->Sort();
533 }
534
535 //________________________________________________________________________
536 Bool_t AliEmcalJetTask::GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const
537 {
538   // Get the leading jets.
539
540   static Float_t pt[9999] = {0};
541
542   const Int_t n = (Int_t)array.size();
543
544   if (n < 1)
545     return kFALSE;
546   
547   for (Int_t i = 0; i < n; i++) 
548     pt[i] = array[i].perp();
549
550   TMath::Sort(n, pt, indexes);
551
552   return kTRUE;
553 }
554
555 //________________________________________________________________________
556 Bool_t AliEmcalJetTask::DoInit()
557 {
558   // Init. Return true if successful.
559
560   if (fTrackEfficiency < 1.) {
561     if (gRandom) delete gRandom;
562     gRandom = new TRandom3(0);
563   }
564
565   // get input collections
566   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
567
568   // get the event
569   fEvent = InputEvent();
570   if (!fEvent) {
571     AliError(Form("%s: Could not retrieve event! Returning", GetName()));
572     return 0;
573   }
574
575   // add jets to event if not yet there
576   if (!(fEvent->FindListObject(fJetsName)))
577     fEvent->AddObject(fJets);
578   else {
579     AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
580     return 0;
581   }
582
583   if (fTracksName == "Tracks")
584     am->LoadBranch("Tracks");
585   if (!fTracks && !fTracksName.IsNull()) {
586     fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
587     if (!fTracks) {
588       AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
589       return 0;
590     }
591   }
592   if (fTracks) {
593     TClass cls(fTracks->GetClass()->GetName());
594     if (cls.InheritsFrom("AliMCParticle") || cls.InheritsFrom("AliAODMCParticle"))
595       fIsMcPart = 1;
596   }
597   
598   if (fCaloName == "CaloClusters")
599     am->LoadBranch("CaloClusters");
600   if (!fClus && !fCaloName.IsNull()) {
601     fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
602     if (!fClus) {
603       AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
604       return 0;
605     }
606   }
607   if (fClus) {
608     TClass cls(fClus->GetClass()->GetName());
609     if (cls.InheritsFrom("AliEmcalParticle"))
610       fIsEmcPart = 1;
611   }
612
613   return 1;
614 }