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