updates to AliEmcalJet, AliEmcalJetTask and AliFJWrapper to include subtraction of...
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliEmcalJetTask.cxx
1 //
2 // Emcal jet finder task.
3 //
4 // Authors: C.Loizides, S.Aiola, M. Verweij
5
6 #include <vector>
7 #include "AliEmcalJetTask.h"
8
9 #include <TChain.h>
10 #include <TClonesArray.h>
11 #include <TList.h>
12 #include <TLorentzVector.h>
13 #include <TMath.h>
14 #include <TRandom3.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 #include "AliRhoParameter.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   fJetsSubName(""),
41   fJetType(kNone),
42   fConstSel(0),
43   fMCConstSel(0),
44   fMarkConst(kFALSE),
45   fRadius(0.4),
46   fMinJetTrackPt(0.15),
47   fMinJetClusPt(0.15),
48   fPhiMin(0),
49   fPhiMax(TMath::TwoPi()),
50   fEtaMin(-0.9),
51   fEtaMax(+0.9),
52   fMinJetArea(0.005),
53   fMinJetPt(1.0),
54   fJetPhiMin(-10),
55   fJetPhiMax(+10),
56   fJetEtaMin(-1),
57   fJetEtaMax(+1),
58   fGhostArea(0.005),
59   fMinMCLabel(0),
60   fRecombScheme(fastjet::pt_scheme),
61   fTrackEfficiency(1.),
62   fMCFlag(0),
63   fGeneratorIndex(-1),
64   fIsInit(0),
65   fLocked(0),
66   fIsPSelSet(0),
67   fIsMcPart(0),
68   fIsEmcPart(0),
69   fLegacyMode(kFALSE),
70   fCodeDebug(kFALSE),
71   fDoGenericSubtractionJetMass(kFALSE),
72   fDoGenericSubtractionGR(kFALSE),
73   fDoConstituentSubtraction(kFALSE),
74   fUseExternalBkg(kFALSE),
75   fRhoName(""),
76   fRhomName(""),
77   fRho(0),
78   fRhom(0),
79   fRMax(0.4),
80   fDRStep(0.04),
81   fPtMinGR(40.),
82   fJets(0),
83   fJetsSub(0),
84   fEvent(0),
85   fTracks(0),
86   fClus(0),
87   fRhoParam(0),
88   fRhomParam(0)
89 {
90   // Default constructor.
91 }
92
93 //________________________________________________________________________
94 AliEmcalJetTask::AliEmcalJetTask(const char *name) : 
95   AliAnalysisTaskSE(name),
96   fTracksName("Tracks"),
97   fCaloName("CaloClusters"),
98   fJetsName("Jets"),
99   fJetsSubName(""),
100   fJetType(kAKT|kFullJet|kRX1Jet),
101   fConstSel(0),
102   fMCConstSel(0),
103   fMarkConst(kFALSE),
104   fRadius(0.4),
105   fMinJetTrackPt(0.15),
106   fMinJetClusPt(0.15),
107   fPhiMin(0),
108   fPhiMax(TMath::TwoPi()),
109   fEtaMin(-0.9),
110   fEtaMax(+0.9),
111   fMinJetArea(0.001),
112   fMinJetPt(1.0),
113   fJetPhiMin(-10),
114   fJetPhiMax(+10),
115   fJetEtaMin(-1),
116   fJetEtaMax(+1),
117   fGhostArea(0.005),
118   fMinMCLabel(0),
119   fRecombScheme(fastjet::pt_scheme),
120   fTrackEfficiency(1.),
121   fMCFlag(0),
122   fGeneratorIndex(-1),
123   fIsInit(0),
124   fLocked(0),
125   fIsPSelSet(0),
126   fIsMcPart(0),
127   fIsEmcPart(0),
128   fLegacyMode(kFALSE),
129   fCodeDebug(kFALSE),
130   fDoGenericSubtractionJetMass(kFALSE),
131   fDoGenericSubtractionGR(kFALSE),
132   fDoConstituentSubtraction(kFALSE),
133   fUseExternalBkg(kFALSE),
134   fRhoName(""),
135   fRhomName(""),
136   fRho(0),
137   fRhom(0),
138   fRMax(0.4),
139   fDRStep(0.04),
140   fPtMinGR(40.),
141   fJets(0),
142   fJetsSub(0),
143   fEvent(0),
144   fTracks(0),
145   fClus(0),
146   fRhoParam(0),
147   fRhomParam(0)
148 {
149   // Standard constructor.
150
151   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
152 }
153
154 //________________________________________________________________________
155 AliEmcalJetTask::~AliEmcalJetTask()
156 {
157   // Destructor
158 }
159
160 //________________________________________________________________________
161 void AliEmcalJetTask::UserCreateOutputObjects()
162 {
163   // Create user objects.
164
165   fJets = new TClonesArray("AliEmcalJet");
166   fJets->SetName(fJetsName);
167
168   if(!fJetsSubName.IsNull()) {
169     fJetsSub = new TClonesArray("AliEmcalJet");
170     fJetsSub->SetName(fJetsSubName);
171   }
172 }
173
174 //________________________________________________________________________
175 void AliEmcalJetTask::UserExec(Option_t *) 
176 {
177   // Main loop, called for each event.
178   if (!fIsInit) {
179     if (!DoInit())
180       return;
181     fIsInit = kTRUE;
182   }
183
184   // clear the jet array (normally a null operation)
185   fJets->Delete();
186   if(fJetsSub) fJetsSub->Delete();
187
188   FindJets();
189 }
190
191 //________________________________________________________________________
192 void AliEmcalJetTask::Terminate(Option_t *) 
193 {
194   // Called once at the end of the analysis.
195 }
196
197 //________________________________________________________________________
198 void AliEmcalJetTask::FindJets()
199 {
200   // Find jets.
201   if (!fTracks && !fClus){
202     cout << "WARNING NO TRACKS OR CLUSTERS:"  <<endl;
203     return;
204   }
205
206  if (fRhoParam)
207     fRho = fRhoParam->GetVal();
208   if (fRhomParam)
209     fRhom = fRhomParam->GetVal();
210
211   TString name("kt");
212   fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
213   if ((fJetType & kAKT) != 0) {
214     name  = "antikt";
215     jalgo = fastjet::antikt_algorithm;
216     AliDebug(1,"Using AKT algorithm");
217   }
218   else {
219     AliDebug(1,"Using KT algorithm");
220   }
221
222   if ((fJetType & kR020Jet) != 0)
223     fRadius = 0.2;
224   else if ((fJetType & kR030Jet) != 0)
225     fRadius = 0.3; 
226   else if ((fJetType & kR040Jet) != 0)
227     fRadius = 0.4;
228
229   // setup fj wrapper
230   AliFJWrapper fjw(name, name);
231   fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
232   fjw.SetGhostArea(fGhostArea);
233   fjw.SetR(fRadius);
234   fjw.SetAlgorithm(jalgo);
235   fjw.SetRecombScheme(static_cast<fastjet::RecombinationScheme>(fRecombScheme)); 
236   fjw.SetMaxRap(fEtaMax);
237   fjw.Clear();
238
239   // get primary vertex
240   Double_t vertex[3] = {0, 0, 0};
241   if(fEvent->GetPrimaryVertex()) fEvent->GetPrimaryVertex()->GetXYZ(vertex);
242
243   AliDebug(2,Form("Jet type = %d", fJetType));
244
245   if ((fIsMcPart || ((fJetType & kFullJet) != 0) || ((fJetType & kChargedJet) != 0)) && fTracks) {
246     const Int_t Ntracks = fTracks->GetEntries();
247     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
248       AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(iTracks));
249       if (!t)
250         continue;
251       if (fIsMcPart) {
252         if (((fJetType & kChargedJet) != 0) && (t->Charge() == 0)) {
253           AliDebug(2,Form("Skipping track %d because it is neutral.", iTracks));
254           continue;
255         }
256         if (((fJetType & kNeutralJet) != 0) && (t->Charge() != 0)) {
257           AliDebug(2,Form("Skipping track %d because it is charged.", iTracks));
258           continue;
259         }
260       }
261       if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) {
262         if (t->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
263           AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask)));
264           continue;
265         }
266         else {
267           AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask)));
268         }
269       }
270       else {
271         if (t->TestBits(fConstSel) != (Int_t)fConstSel) {
272           AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));
273           continue;
274         }
275         else {
276           AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));           
277         }
278       }
279       if ((t->GetFlag() & fMCFlag) != fMCFlag) { 
280         AliDebug(2,Form("Skipping track %d because it does not match the MC flags", iTracks));
281         continue;
282       }
283       if (fGeneratorIndex >= 0 && t->GetGeneratorIndex() != fGeneratorIndex) { 
284         AliDebug(2,Form("Skipping track %d because it does not match the MC generator index", iTracks));
285         continue;
286       }
287       if (t->Pt() < fMinJetTrackPt) 
288         continue;
289       Double_t eta = t->Eta();
290       Double_t phi = t->Phi();
291       if ((eta<fEtaMin) || (eta>fEtaMax) ||
292           (phi<fPhiMin) || (phi>fPhiMax))
293         continue;
294
295       // artificial inefficiency
296       if (fTrackEfficiency < 1.) {
297         Double_t rnd = gRandom->Rndm();
298         if (fTrackEfficiency < rnd) {
299           AliDebug(2,Form("Track %d rejected due to artificial tracking inefficiency", iTracks));
300           continue;
301         }
302       }
303
304       // offset of 100 for consistency with cluster ids
305       AliDebug(2,Form("Track %d accepted (label = %d, pt = %f)", iTracks, t->GetLabel(), t->Pt()));
306       fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->E(), iTracks + 100);
307     }
308   }
309
310   if ((((fJetType & kFullJet) != 0) || ((fJetType & kNeutralJet) != 0)) && fClus) {
311     const Int_t Nclus = fClus->GetEntries();
312     for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
313       AliVCluster *c = 0;
314       Double_t cEta=0,cPhi=0,cPt=0;
315       Double_t cPx=0,cPy=0,cPz=0;
316       if (fIsEmcPart) {
317         AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(iClus));
318         if (!ep)
319           continue;
320
321         c = ep->GetCluster();
322         if (!c)
323           continue;
324
325         if (c->GetLabel() > fMinMCLabel) {
326           if (ep->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
327             AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
328             continue;
329           }
330           else {
331             AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
332           }
333         }
334         else {
335           if (ep->TestBits(fConstSel) != (Int_t)fConstSel) {
336             AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));
337             continue;
338           }
339           else {
340             AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));        
341           }
342         }
343
344         cEta = ep->Eta();
345         cPhi = ep->Phi();
346         cPt  = ep->Pt();
347         cPx  = ep->Px();
348         cPy  = ep->Py();
349         cPz  = ep->Pz();
350       } else {
351         c = static_cast<AliVCluster*>(fClus->At(iClus));
352         if (!c)
353           continue;
354
355         if (c->GetLabel() > fMinMCLabel) {
356           if (c->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
357             AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
358             continue;
359           }
360           else {
361             AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
362           }
363         }
364         else {
365           if (c->TestBits(fConstSel) != (Int_t)fConstSel) {
366             AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));
367             continue;
368           }
369           else {
370             AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));         
371           }
372         }
373
374         TLorentzVector nP;
375         c->GetMomentum(nP, vertex);
376         cEta = nP.Eta();
377         cPhi = nP.Phi();
378         cPt  = nP.Pt();
379         cPx  = nP.Px();
380         cPy  = nP.Py();
381         cPz  = nP.Pz();
382       }
383       if (!c->IsEMCAL())
384         continue;
385       if (cPt < fMinJetClusPt) 
386         continue;
387       if ((cEta<fEtaMin) || (cEta>fEtaMax) ||
388           (cPhi<fPhiMin) || (cPhi>fPhiMax))
389         continue;
390       // offset of 100 to skip ghost particles uid = -1
391       AliDebug(2,Form("Cluster %d accepted (label = %d)", iClus, c->GetLabel()));
392       fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz), -iClus - 100);
393     }
394   }
395
396   // setting legacy mode
397   if (fLegacyMode) { 
398     fjw.SetLegacyMode(kTRUE);
399   }
400
401   // run jet finder
402   fjw.Run();
403
404   //run generic subtractor
405   if(fDoGenericSubtractionJetMass) {
406     fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
407     fjw.DoGenericSubtractionJetMass();
408   }
409
410   //run constituent subtractor
411   if(fDoConstituentSubtraction) {
412     fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
413     fjw.DoConstituentSubtraction();
414   }
415
416   // get geometry
417   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
418   if (!geom) {
419     AliFatal(Form("%s: Can not create geometry", GetName()));
420     return;
421   }
422
423   // loop over fastjet jets
424   std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
425   // sort jets according to jet pt
426   static Int_t indexes[9999] = {-1};
427   GetSortedArray(indexes, jets_incl);
428
429   AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
430   for (UInt_t ijet=0, jetCount=0; ijet<jets_incl.size(); ++ijet) {
431     Int_t ij = indexes[ijet];
432     AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fjw.GetJetArea(ij)));
433
434     if (jets_incl[ij].perp()<fMinJetPt) 
435       continue;
436     if (fjw.GetJetArea(ij)<fMinJetArea)
437       continue;
438     if ((jets_incl[ij].eta()<fJetEtaMin) || (jets_incl[ij].eta()>fJetEtaMax) ||
439         (jets_incl[ij].phi()<fJetPhiMin) || (jets_incl[ij].phi()>fJetPhiMax))
440       continue;
441
442     AliEmcalJet *jet = new ((*fJets)[jetCount]) 
443       AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
444     jet->SetLabel(ij);
445
446     //do generic subtraction if requested
447 #ifdef FASTJET_VERSION
448     if(fDoGenericSubtractionJetMass) {
449       std::vector<fastjet::contrib::GenericSubtractorInfo> jetMassInfo = fjw.GetGenSubtractorInfoJetMass();
450       Int_t n = (Int_t)jetMassInfo.size();
451       if(n>ij && n>0) {
452         jet->SetFirstDerivative(jetMassInfo[ij].first_derivative());
453         jet->SetSecondDerivative(jetMassInfo[ij].second_derivative());
454         jet->SetFirstOrderSubtracted(jetMassInfo[ij].first_order_subtracted());
455         jet->SetSecondOrderSubtracted(jetMassInfo[ij].second_order_subtracted());
456       }
457     }
458     //here do generic subtraction for angular structure function
459     Double_t ptcorr = jets_incl[ij].perp()-fjw.GetJetArea(ij)*fRho;
460     fRMax = fRadius+0.2;
461     if(fDoGenericSubtractionGR && ptcorr>fPtMinGR) {
462       fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
463       fjw.SetRMaxAndStep(fRMax,fDRStep);
464       fjw.DoGenericSubtractionGR(ij);
465       std::vector<double> num = fjw.GetGRNumerator();
466       std::vector<double> den = fjw.GetGRDenominator();
467       std::vector<double> nums = fjw.GetGRNumeratorSub();
468       std::vector<double> dens = fjw.GetGRDenominatorSub();
469       //pass this to AliEmcalJet
470       jet->SetGRNumSize(num.size());
471       jet->SetGRDenSize(den.size());
472       jet->SetGRNumSubSize(nums.size());
473       jet->SetGRDenSubSize(dens.size());
474       Int_t nsize = (Int_t)num.size();
475       for(Int_t g = 0; g<nsize; ++g) {
476         jet->AddGRNumAt(num[g],g);
477         jet->AddGRNumSubAt(nums[g],g);
478       }
479       Int_t dsize = (Int_t)den.size();
480       for(Int_t g = 0; g<dsize; ++g) {
481         jet->AddGRDenAt(den[g],g);
482         jet->AddGRDenSubAt(dens[g],g);
483       }
484     }
485 #endif
486
487     // loop over constituents
488     std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
489     jet->SetNumberOfTracks(constituents.size());
490     jet->SetNumberOfClusters(constituents.size());
491
492     Int_t nt            = 0;
493     Int_t nc            = 0;
494     Double_t neutralE   = 0;
495     Double_t maxCh      = 0;
496     Double_t maxNe      = 0;
497     Int_t gall          = 0;
498     Int_t gemc          = 0;
499     Int_t cemc          = 0;
500     Int_t ncharged      = 0;
501     Int_t nneutral      = 0;
502     Double_t mcpt       = 0;
503     Double_t emcpt      = 0;
504
505     for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
506       Int_t uid = constituents[ic].user_index();
507       AliDebug(3,Form("Processing constituent %d", uid));
508       if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
509         ++gall;
510         Double_t gphi = constituents[ic].phi();
511         if (gphi<0) 
512           gphi += TMath::TwoPi();
513         gphi *= TMath::RadToDeg();
514         Double_t geta = constituents[ic].eta();
515         if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
516             (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
517           ++gemc; 
518       } else if ((uid > 0) && fTracks) { // track constituent
519         Int_t tid = uid - 100;
520         AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
521         if (!t) {
522           AliError(Form("Could not find track %d",tid));
523           continue;
524         }
525         if (jetCount < fMarkConst) {
526           AliDebug(2,Form("Marking track %d with bit map %d", tid, fJetType));
527           t->SetBit(fJetType);
528         }
529         Double_t cEta = t->Eta();
530         Double_t cPhi = t->Phi();
531         Double_t cPt  = t->Pt();
532         Double_t cP   = t->P();
533         if (t->Charge() == 0) {
534           neutralE += cP;
535           ++nneutral;
536           if (cPt > maxNe)
537             maxNe = cPt;
538         } else {
539           ++ncharged;
540           if (cPt > maxCh)
541             maxCh = cPt;
542         }
543         if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) // check if MC particle
544           mcpt += cPt;
545
546         if (cPhi<0) 
547           cPhi += TMath::TwoPi();
548         cPhi *= TMath::RadToDeg();
549         if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
550             (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
551           emcpt += cPt;
552           ++cemc;
553         }
554
555         jet->AddTrackAt(tid, nt);
556         ++nt;
557       } else if (fClus) { // cluster constituent
558         Int_t cid = -uid - 100;
559         AliVCluster *c = 0;
560         Double_t cEta=0,cPhi=0,cPt=0,cP=0;
561         if (fIsEmcPart) {
562           AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
563           if (!ep)
564             continue;
565           c = ep->GetCluster();
566           if (!c)
567             continue;
568           if (jetCount < fMarkConst)
569             ep->SetBit(fJetType);
570           cEta = ep->Eta();
571           cPhi = ep->Phi();
572           cPt  = ep->Pt();
573           cP   = ep->P();
574         } else {
575           c = static_cast<AliVCluster*>(fClus->At(cid));
576           if (!c)
577             continue;
578           if (jetCount < fMarkConst)
579             c->SetBit(fJetType);
580           TLorentzVector nP;
581           c->GetMomentum(nP, vertex);
582           cEta = nP.Eta();
583           cPhi = nP.Phi();
584           cPt  = nP.Pt();
585           cP   = nP.P();
586         }
587
588         neutralE += cP;
589         if (cPt > maxNe)
590           maxNe = cPt;
591
592         if (c->GetLabel() > fMinMCLabel) // MC particle
593           mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
594
595         if (cPhi<0) 
596           cPhi += TMath::TwoPi();
597         cPhi *= TMath::RadToDeg();
598         if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
599             (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
600           emcpt += cPt;
601           ++cemc;
602         }
603
604         jet->AddClusterAt(cid, nc);
605         ++nc;
606         ++nneutral;
607       } else {
608         AliError(Form("%s: No logical way to end up here.", GetName()));
609         continue;
610       }
611     } /* end of constituent loop */
612
613     jet->SetNumberOfTracks(nt);
614     jet->SetNumberOfClusters(nc);
615     jet->SortConstituents();
616     jet->SetMaxNeutralPt(maxNe);
617     jet->SetMaxChargedPt(maxCh);
618     jet->SetNEF(neutralE / jet->E());
619     fastjet::PseudoJet area(fjw.GetJetAreaVector(ij));
620     jet->SetArea(area.perp());
621     jet->SetAreaEta(area.eta());
622     jet->SetAreaPhi(area.phi());
623     jet->SetNumberOfCharged(ncharged);
624     jet->SetNumberOfNeutrals(nneutral);
625     jet->SetMCPt(mcpt);
626     jet->SetNEmc(cemc);
627     jet->SetPtEmc(emcpt);
628
629     if (gall > 0)
630       jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
631     else 
632       jet->SetAreaEmc(-1);
633     if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) && 
634         (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
635         (jet->Eta() > geom->GetArm1EtaMin()) && 
636         (jet->Eta() < geom->GetArm1EtaMax()))
637       jet->SetAxisInEmcal(kTRUE);
638
639     AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), (Int_t)constituents.size()));
640     jetCount++;
641   }
642   //fJets->Sort();
643
644   //run constituent subtractor if requested
645   if(fDoConstituentSubtraction) {
646 #ifdef FASTJET_VERSION
647     if(!fJetsSub) AliWarning(Form("No jet branch to write to for subtracted jets. fJetsSubName: %s",fJetsSubName.Data()));
648     else {
649       std::vector<fastjet::PseudoJet> jets_sub;
650       jets_sub = fjw.GetConstituentSubtrJets();
651       AliDebug(1,Form("%d constituent subtracted jets found", (Int_t)jets_sub.size()));
652       for (UInt_t ijet=0, jetCount=0; ijet<jets_sub.size(); ++ijet) {
653         //Only storing 4-vector and jet area of unsubtracted jet        
654         AliEmcalJet *jet_sub = new ((*fJetsSub)[jetCount]) 
655           AliEmcalJet(jets_sub[ijet].perp(), jets_sub[ijet].eta(), jets_sub[ijet].phi(), jets_sub[ijet].m());
656         jet_sub->SetLabel(ijet);
657         fastjet::PseudoJet area(fjw.GetJetAreaVector(ijet));
658         jet_sub->SetArea(area.perp());
659         jet_sub->SetAreaEta(area.eta());
660         jet_sub->SetAreaPhi(area.phi());
661         jet_sub->SetAreaEmc(area.perp());
662         jetCount++;
663       }
664     }
665 #endif
666   } //constituent subtraction
667 }
668
669 //________________________________________________________________________
670 Bool_t AliEmcalJetTask::GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const
671 {
672   // Get the leading jets.
673
674   static Float_t pt[9999] = {0};
675
676   const Int_t n = (Int_t)array.size();
677
678   if (n < 1)
679     return kFALSE;
680   
681   for (Int_t i = 0; i < n; i++) 
682     pt[i] = array[i].perp();
683
684   TMath::Sort(n, pt, indexes);
685
686   return kTRUE;
687 }
688
689 //________________________________________________________________________
690 Bool_t AliEmcalJetTask::DoInit()
691 {
692   // Init. Return true if successful.
693
694   if (fTrackEfficiency < 1.) {
695     if (gRandom) delete gRandom;
696     gRandom = new TRandom3(0);
697   }
698
699   // get input collections
700   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
701
702   // get the event
703   fEvent = InputEvent();
704   if (!fEvent) {
705     AliError(Form("%s: Could not retrieve event! Returning", GetName()));
706     return 0;
707   }
708
709   // add jets to event if not yet there
710   if (!(fEvent->FindListObject(fJetsName)))
711     fEvent->AddObject(fJets);
712   else {
713     AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
714     return 0;
715   }
716
717   if (!(fEvent->FindListObject(fJetsSubName)) && fJetsSub)
718     fEvent->AddObject(fJetsSub);
719
720   if (fTracksName == "Tracks")
721     am->LoadBranch("Tracks");
722   if (!fTracks && !fTracksName.IsNull()) {
723     fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
724     if (!fTracks) {
725       AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
726       return 0;
727     }
728   }
729   if (fTracks) {
730     TClass cls(fTracks->GetClass()->GetName());
731     if (cls.InheritsFrom("AliMCParticle") || cls.InheritsFrom("AliAODMCParticle"))
732       fIsMcPart = 1;
733   }
734   
735   if (fCaloName == "CaloClusters")
736     am->LoadBranch("CaloClusters");
737   if (!fClus && !fCaloName.IsNull()) {
738     fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
739     if (!fClus) {
740       AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
741       return 0;
742     }
743   }
744   if (fClus) {
745     TClass cls(fClus->GetClass()->GetName());
746     if (cls.InheritsFrom("AliEmcalParticle"))
747       fIsEmcPart = 1;
748   }
749
750   if (!fRhoName.IsNull() && !fRhoParam) { // get rho from the event
751     fRhoParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
752     if (!fRhoParam) {
753       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
754       return 0;
755     }
756   }
757   if (!fRhomName.IsNull() && !fRhomParam) { // get rhom from the event
758     fRhomParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhomName));
759     if (!fRhomParam) {
760       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhomName.Data()));
761       return 0;
762     }
763   }
764   
765   return 1;
766 }