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